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Existing  routing  models  for  hydrologic  systems  have  either  1)  a 
complex,  nonlinear  structure  not  amenable  to  second  order  (i.e.,  error) 
analysis,  requiring  considerable  computational  resources  for  solution; 
or  2)  a  simple,  linear  structure  that  is  strictly  valid  over  only  a 
small  flow  range  about  which  the  model  is  linearized. 

The  new  perspective  for  routing  surface  runoff  or  channel  flow 
uses  any  appropriate  model,  called  the  primary  model,  to  simulate  the 
observed  physics  for  the  catchment  or  channel  of  interest,  and  then 
emulates  the  mathematical  behavior  of  the  primary  model  with  a 
computationally  efficient  nonlinear  state  space  model  that  has 
structural  properties  amenable  to  second  order  analysis.   The  overall 
strategy  matches  impulse  response  properties  of  the  state  space  model  to 
those  of  the  primary  model  over  a  wide  range  of  flow  levels.   In 
general,  impulse  response  properties  vary  with  flow  level. 
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Two  important  response  properties  are  time  lag  and  attenuation. 
Nonlinearities  in  time  lag  occur  because  water  flows  faster  as  it  flows 
deeper.   Nonlinearities  in  attenuation  occur  because  deep  waves  disperse 
more  rapidly  than  shallow  waves.   These  two  properties  can  be  expressed 
as  linear  or  piecewise  linear  functions  of  discharge  in  log-space,  thus 
relating  them  directly  to  the  parameters  of  a  nonlinear  state  space 
model. 

The  limitations  of  this  approach  to  routing  were  examined  by 
comparing  results  from  the  nonlinear  state  space  model  with  a  numerical 
solution  of  the  Saint-Venant  equations.  No  significant  differences  were 
observed.   The  state  space  model  required  two  to  three  orders  of 
magnitude  less  CPU  time  than  numerical  solution  of  the  Saint-Venant 
equations . 

A.  generalized,  computationally  efficient  state  space  routing  model 
that  can  emulate  the  nonlinear  behavior  of  any  other  routing  model  for 
downstream  waves  has  been  developed.   Further  studies  are  required  to 
determine  the  full  usefulness  of  this  modelling  approach. 
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CHAPTER  I 
A  NEW  PERSPECTIVE  OF  NONLINEAR  ROUTING 


Statement  of  the  Problem 

Existing  technology  for  routing  unsteady  flow  in  open  channels  and 
for  routing  surface  runoff  in  a  catchment  has  not  proved  adequate  for 
some  important  hydrologic  applications.   Current  routing  methods  are 
either  inefficient  when  numerous  simulations  are  required,  as  in  flood 
frequency  analysis;  inadequate  when  the  forecast  requires  an  estimate  of 
uncertainty  in  past  or  future  flow  conditions,  as  in  water  supply 
forecasting;  or  unsatisfactory  when  the  information  content  of  current 
data  is  important,  as  in  real-time  river  forecasting.   The  above 
applications  require 

1.  an  objective  estimate  of  uncertainty, 

2.  repetitive  computations,  and/or 

3.  consideration  of  second  order  (i.e.,  error)  properties  of 
models  and  data  as  well  as  the  first  order  (i.e.,  mean) 
trajectory  in  time  and  space  of  the  flow  conditions; 

as  well  as  requiring  models  or  computational  algorithms  to  account  for 
significant  nonlinearities  of  the  channel  or  catchment  response. 

For  example,  the  National  Weather  Service  (NWS),  the  Bureau  of 
Reclamation,  and  the  Federal  Emergency  Management  Authority  typically 
study  the  thousands  of  dams  throughout  the  country  by  simulating  various 
dam  failure  scenarios.   Numerical  solution  of  the  full  Saint-Venant 
equations  can  impose  a  significant  computational  burden  in  such 


contingency  dam  break  simulations  for  which  repetitive  calculations  are 
required.   Simpler  models,  while  computationally  efficient,  do  not 
properly  account  for  the  nonlinear  behavior  important  to  these  analyses. 

The  NWS  studies  of  various  possible  dam  failures  are  designed  to 
provide  timely  and  accurate  forecasts  of  expected  flow  levels  to 
communities  downstream  of  the  dam  in  the  event  of  a  dam  failure.   When 
the  threat  of  a  dam  failure  is  imminent,  the  NWS  Forecast  Office  can 
provide  worst  case  (i.e.,  instantaneous  failure)  and  alternative 
projected  flow  conditions  to  the  media  and  emergency  authorities. 
Several  computer  simulations  involving  flow  routing  downstream  from  each 
dam  studied  are  required  to  produce  that  information. 

For  daily,  real-time  river  forecasting,  as  performed  at  NWS  River 
Forecast  Centers,  timely  processing  of  all  current  observations  is 
critical.   A  rapid  mechanism  is  needed  to  update  the  conditions 
projected  by  the  computer  model  to  account  for  the  most  recently 
observed  data.   A  state  space  model  of  the  full  Saint-Venant  equations, 
which  would  be  amenable  to  updating  with  modern  estimation  theory, 
imposes  a  large  computational  burden.   Simpler  routing  models  are  not 
hydraulically  adequate  over  the  range  of  inflow  transients  encountered 
in  flood  forecasting.   Although  some  simple  models  could  be  formulated 
for  use  with  updating  algorithms,  their  forecasting  capability  is 
diminished. by  the  limited  flow  range  over  which  they  properly  model  the 
nonlinear  nature  of  water  movement. 

Flood  frequency  studies,  such  as  those  performed  during  the 
hydraulic  design  of  channel  structures,  often  require  computer 
simulations  for  many  events  over  a  wide  range  of  flow  conditions. 
Again,  numerical  solution  of  the  full  Saint-Venant  equations  is  costly 


and  the  simpler,  linear  models  do  not  properly  simulate  for  large 
variations  of  inflow. 

Statement  of  the  New  Perspective 

The  new  perspective  of  routing  in  hydrology  is  to  let  any  model, 
called  the  primary  model,  simulate  the  physics  of  a  channel  or 
catchment,  and  to  emulate  the  behavior  of  the  primary  model  with  a  model 
structured  to  meet  the  applications  described  above.   This  process  is 
depicted  in  Figure  1.1.   To  account  for  the  modelling  of  second  order 
properties  (i.e.,  variance  and  covariance  measurements  of  uncertainty) 
and  to  simplify  computation,  a  state  space  form  has  been  chosen  as  the 
model  to  emulate  primary  model  behavior.   The  state  space  model 
structure  is 

X    =F«X+G„-U,  (I'D 

-t+1   -X,t  -t   -X,t  -t+1 

where  X^        =  the  system  state  vector, 

Fv  f  =  the  state  transition  matrix, 

A.,  t 

G^  t  =  the  input  coefficient  matrix,  and 
Ur+1  =  the  vector  of  system  inputs. 

The  structure  of  a  state  space  model  is  possibly  the  only  one  that 

meets  the  previously  described  objective  of  simplicity  and  has  the 

potential  for  simulating  second  order  properties.   The  state  space  model 

can  be  used  with  modern  estimation  theory  in  a  Kalman  filtering 

algorithm  to  propagate  and  update  both  mean  and  variance  properties  of 

the  channel  or  catchment.   The  primary  model  simulates  the  physics  of 

the  hydrologic  system  while  the  state  space  model,  because  of  its 

desired  computational  and  structural  properties,  emulates  the 

mathematical  behavior  of  the  primary  model. 
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Figure  1.1    An  Overview  of  the  New  Perspective 
of  Nonlinear  Routing  in  Hydrology 


The  essence  of  this  new  perspective  of  routing  is  its  two-step 
approach  to  the  problems  encountered  with  existing  routing  technology. 
First,  the  physics  of  the  system  are  modelled  with  a  primary  model,  then 
the  mathematics  of  the  primary  model  are  emulated  by  a  state  space 
model.   The  state  space  model  is  a  model  of  the  primary  model. 

Development  of  the  New  Perspective 
Overview 


Several  steps  are  required  to  apply  this  new  perspective  to  a 
routing  problem  of  interest.   The  steps  are 

1.  select  a  primary  model, 

2.  calibrate  the  primary  model  to  the  hydrologic  system, 

3.  determine  the  appropriate  structure  for  the  state  space 
model, 

4.  calibrate  the  state  space  model  to  the  primary  model,  and 

5.  emulate  the  primary  model  behavior  with  the  state  space 
model. 

Selection  of  a  Primary  Model 

The  primary  model  simulates  the  physics  of  a  channel  or 
catchment.   Because,  in  general,  the  forces  governing  the  motion  of 
water  interact  in  a  complicated  manner,  a  complex  model  typically  will 
be  most  appropriate  to  describe  the  physical  properties  of  the 
hydrologic  system.   However,  the  primary  model  does  not  have  to  be 
complex.   Any  simple  model  can  be  chosen  to  simulate  the  system  physics 
based  on  data  or  computational  limitations,  or  modeller  preference.   A 


summary  of  some  available  channel  routing  models  is  presented  in 
Appendix  A. 

Calibration  of  the  Primary  Model  to  the  Hydrologic  System 

The  effort  required  to  calibrate  the  primary  model  will  vary 
greatly  depending  upon  the  complexity  of  the  primary  model  chosen.   If  a 
physically  detailed  hydraulic  model  is  selected,  channel  cross-section 
data  must  be  obtained  and  reduced  to  the  form  required  by  the  model. 
This  requires  a  complicated  process  of  adjusting  friction  factors, 
determining  active  versus  inactive  flow  area,  etc.   For  simpler  models, 
the  calibration  process  usually  requires  determining  a  few  parameters 
based  on  the  comparison  of  simulated  and  observed  data  for  some  past 
inflows.   Parameter  selection  for  simple  models  can  often  be  automated, 
given  inflow  and  outflow  records.   Even  for  some  complex  hydraulic 
models,  final  parameter  adjustment  can  be  performed  automatically  (Fread 
and  Smith  1978). 

Structure  of  the  State  Space  Model 

The  mathematical  properties  of  the  primary  model  can  be 
represented  in  numerous  ways,  given  the  general  structure  of  equation 
1.1.   The  system  states  (i.e.,  the  elements  of  the  vector  X_)  can 
represent  many  possible  quantities.   The  states  of  the  surrogate  model 
are  determined  by  the  important  relationships  in  the  primary  model. 
State  space  models  of  hydrologic  systems  can  be  constructed  with  ^ 
physically  based  states  such  as  the  quantity  of  water  in  a  reservoir 
(Szollosi-Nagy  1981)  or  with  states  which  are  purely  mathematical 


artifacts  that  reproduce  observed  system  behavior  (Goldstein  and 
Larimore  1980).   In  fact,  the  full  Saint-Venant  equations  can  be 
expressed  in  state  space  form.   At  each  time  step,  values  of  flow  and 
cross-sectional  area  at  numerous  points  along  the  channel  represent  the 
state  of  the  system.   The  state  vector,  Xj. ,  could  be  constructed  from 
these  flow  and  area  values.   However,  since  the  number  of  points  defined 
along  a  channel  typically  is  large,  the  state  vector  would  be 
correspondingly  large.   For  example,  some  NWS  applications  of  the 
Dynamic  Wave  Operational  Program  (DWOPER),  which  numerically  solves  the 
Saint-Venant  equations  (Fread  and  Smith  1978)  have  25  and  45  points 
defined  along  the  channels,  respectively,  for  the  Columbia-Willamette 
River  system  and  the  Ohio-Mississippi  River  junction.   The  state  vectors 
for  these  systems  would  consist  of  50  and  90  elements,  respectively, 
with  state  transition  matrices  of  2500  and  8100  elements.   Obvious 
dimensionality  problems  arise  in  the  computer  solution  of  equation  1.1 
for  systems  of  these  sizes. 

The  structure  of  the  nonlinear  state  space  model  developed  is 
presented  in  detail  in  Chapter  IV.   The  approach  taken  uses  the 
nonstationary  linear  structure  of  the  state  space  model  as  a  surrogate 
for  the  stationary  nonlinear  behavior  of  the  primary  model. 

Calibration  of  the  State  Space  Model  to  the  Primary  Model 

Problems  encountered  in  certain  routing  applications  can  be  solved 
by  using  the  state  space  model  to  emulate  the  mathematical  behavior  of 
the  primary  model.   Since  the  state  space  model  is  calibrated  to  the 
mathematical  properties  of  the  primary  model,  many  approaches  are 
possible. 


The  approach  taken  here  draws  on  the  general,  powerful  theory  of 
linear  systems.   As  described  in  Chapter  II,  the  nearly  linear  behavior 
of  the  primary  model  for  a  small  flow  range  is  emulated  with  an 
equivalent  linear  system.   This  is  accomplished  by  matching  the  impulse 
responses  of  the  primary  and  state  space  models  for  small  flow  ranges 
about  various  reference  discharges. 

The  inherent  nonlinearities  that  cause  the  behavior  of  the  state 
space  model  to  agree  with  or  to  depart  from  the  primary  model  behavior 
are  variations  with  discharge  of  the  lag  and  dispersive  properties  of 
the  impulse  response  function.   These  properties  vary  substantially  over 
large  changes  of  inflow.   A  simple  functional  relationship  can  be 
derived  that  specifies  the  variation  of  F_  and  ^with  discharge  to 
emulate  properly  the  lag  and  dispersive  properties.   Given  such  a 
functional  relationship,  the  state  space  model  can  emulate  the  full 
nonlinear  and  dynamic  behavior  of  the  primary  model  over  all  flow 
ranges. 

The  State  Space  Model  as  Surrogate  for  the  Primary  Model 

Simple  linear  models  have  been  used  in  place  of  complex  routing 
models  for  many  years.   They  replace  the  complex  model,  but  do  not 
properly  represent  some  of  the  important  properties  of  the  complex 
model.   With  the  new  perspective,  the  primary  model  is  not  replaced. 
Instead  its  mathematical  structure  is  emulated  by  a  state  space  model 
with  computational  and  structural  properties  appropriate  for  the 
solution  of  problems  encountered  in  many  routing  applications.   The 
state  space  model  acts  as  a  surrogate  for  the  primary  model. 


The  state  space  model  typically  has  much  lower  dimensionality  than 
the  primary  model  and  imposes  a  much  smaller  computational  burden. 
Because  the  state  space  model  emulates  the  full  nonlinear  behavior  of 
the  primary  model,  it  can  be  used  in  applications  where  the  primary 
model  is  not  appropriate  for  error  analysis  because  of  the  computational 
burden  or  improper  structure.   This  approach  presents  a  model 
mathematically  equivalent  to  the  primary  model,  that  can  be  used  on 
smaller  computers  than  would  be  possible  with  the  primary  model  itself. 

An  Application  of  the  New  Perspective 

To  explore  the  new  perspective's  applicability  to  routing  in 
hydrology,  the  general  steps  described  above  were  followed  for  a 
specific  primary  model.   As  shown  in  Figure  1.2,  the  application 
selected  was  the  full,  nonlinear  Saint-Venant  equations  as  the  exemplar 
primary  model  to  simulate  downstream  flow  in  a  prismatic  channel.   This 
primary  model  and  physical  system  were  chosen  because 

1.  the  Saint-Venant  equations  are  generally  recognized  as  the 
most  complete  representation  of  unsteady  flow  phenomena, 

2.  an  analytical  solution  for  the  impulse  response  of  the 
Saint-Venant  equations  can  be  derived  for  downstream  wave 
movement  in  prismatic  channels,  and 

3.  numerical  solution  of  the  Saint-Venant  equations  is 
well-studied,  and  many  excellent  solution  algorithms  are 
available  for  comparison  with  the  surrogate  model  (Liggett 
and  Cunge  1975). 

Analysis  of  the  ability  of  a  state  space  model  to  emulate  the 
behavior  of  the  full,  nonlinear  Saint-Venant  equations  should  lend 
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insight  into  the  applicability  of  this  approach  for  other  primary 
models.   Because,  in  general,  no  analytical  solution  for  the  impulse 
response  of  the  primary  model  is  available,  the  case  chosen  allows 
comparison  of  the  true  impulse  response  with  results  obtained 
numerically. 

Organization  of  the  Dissertation 

The  organization  of  this  dissertation  is  schematically  presented 
in  Figure  1.2.   The  new  perspective  of  routing  and  the  reasons  for 
selecting  the  Saint-Venant  equations  as  the  exemplar  primary  model  have 
been  presented.   A  linear  analysis  of  the  Saint-Venant  equations  and 
development  of  a  linear  state  space  model  follow  in  Chapter  II.   In 
Chapter  III,  results  from  the  linear  state  space  model  are  compared  with 
the  full  Saint-Venant  equations  for  small  inflow  transients  about  a 
reference  flow  level.   These  results  indicate  that  the  linear  state 
space  model  is  a  very  adequate  surrogate  for  the  Saint-Venant  equations 
over  small  flow  ranges.   Further  nonlinear  studies  are  motivated  by 
comparison  of  the  linear  state  space  model  with  the  primary  model  for 
large  inflow  hydrographs.   Nonlinear  behavior  of  the  primary  model 
causes  the  comparison  of  results  with  the  linear  state  space  model  to 

deteriorate. 

A  nonlinear  analysis  of  the  Saint-Venant  equations  primary  model 
and  the  development  of  a  nonlinear  state  space  model  are  conducted  in 
Chapter  IV.   Chapter  V  is  devoted  to  the  verification  of  results  for  the 
nonlinear  state  space  model  versus  the  primary  model  for  large  inflow 
hydrographs.   Results  from  the  primary  and  surrogate  models  for  a 
complex  inflow  hydrograph  were  comparable.   However,  the  nonlinear  state 
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space  model  used  two  to  three  orders  of  magnitude  less  CPU  time  than  the 
numerical  solution  of  the  Saint-Venant  equations. 

A  general  procedure  to  emulate  any  primary  model  with  the 
nonlinear  state  space  model  is  presented  in  Chapter  VI.  By  following 
the  steps  here  described,  the  nonlinear  state  space  model  can  be  used  as 
a  surrogate  for  any  channel  routing  or  surface  runoff  model. 

The  relationship  of  the  current  work  to  previous  studies  of 
routing  in  hydrology  is  the  subject  of  Chapter  VII.  A  review  of 
pertinent  work,  by  other  authors,  some  data  to  support  the  power  function 
relationship  between  flow  level  and  time  lag,  and  a  comparison  of  the 
nonlinear  state  space  and  Muskingum-Cunge  models  with  the  full 
Saint-Venant  equations  are  the  major  topics  of  this  chapter. 

Conclusions  on  the  new  perspective  of  routing  in  general,  and  on 
its  application  to  the  exemplar  primary  model  in  particular,  are  stated 
in  Chapter  VIII.   Recommendations  for  future  work,  Including  the 
necessary  tasks  for  the  application  of  estimation  theory  with  the 
nonlinear  state  space  model  and  possible  alternative  approaches  for 
representing  the  nonlinear  model  structure,  form  the  final  portion  of 
this  last  chapter. 


CHAPTER  II 
INSIGHTS  INTO  NONLINEAR  ROUTING  FROM  LINEAR  SYSTEMS  THEORY 


Linear  Systems  Theory 
Introduction 

Linear  systems  theory  provides  many  tools  for  the  analysis  of 
linear  phenomona.   Linear  systems  theory,  strictly  applied,  would  have 
few  uses  because,  in  general,  the  only  truly  linear  systems  are 
mathematical  abstractions.   However,  the  spirit  of  linear  systems 
analysis  lends  itself  to  unraveling  those  aspects  of  nonlinear  systems 
that  behave  in  a  nearly  linear  manner.   The  phenomenon  of  routing  in 
hydrology  is,  in  general,  a  very  nonlinear  process.   Nevertheless, 
important  insights  into  hydrologic  system  behavior  may  be  gained  by 
applying  some  of  the  tools  of  linear  systems  theory. 

The  concepts  of  linear  systems  theory  have  been  applied  to  many 
and  varied  systems  and  are  well  known.   Those  concepts  germane  to  this 
work  are  summarized  below. 

Systems  can  be  thought  of  abstractly  as  the  mechanisms  that 
interrelate  two  objects  (Dooge  1973,  pp.  3-4).   The  operational 
definition  used  in  this  work  is  that  a  system  is  the  process  that 
transforms  input  signals  into  output  signals.   For  the  hydrologic 
processes  analyzed  here,  the  complex  relationships  among  the  forces 
governing  the  motion  of  particles  of  water  are  conceptualized  in  a 
system  as  shown  in  Figure  2.1.   The  value  of  this  view  of  a  system  lies 
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INPUT  SIGNALS 

SYSTEM 

OUTPUT  SIGNALS 

x,U) 

y,«-> 

x2(t) 

• 
• 

a2tt) 

■ 
• 

Figure  2.1  A  System 
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in  the  analytical  tools  it  provides  since  the  mathematical  properties  of 
the  system  determine  which  tools  are  applicable. 

Properties  of  a  System 

The  properties  of  consequence  for  hydrologic  systems  are  whether 
the  system  is 

a.  determinsitic  or  probabilistic, 

b.  distributed  input/distributed  output 

c.  causal, 

d.  stationary,  and 

e.  linear. 

Typically  the  response  of  a  time-varying,  nonstationary,  distributed, 
probabilistic,  nonlinear  hydrologic  system  is  simulated  by  a  time 
invariant,  stationary,  lumped,  deterministic,  linear  model.   And  we 
wonder  why  there  are  discrepancies  between  observations  and  simulated 

outputs! 

A  system  is  determinsitic  if  the  precise  state  of  the  system  (the 
values  of  X  in  equation  1.1)  can  be  foretold.  If  a  random  component  is 
introduced  into  the  variation  of  the  system  states,  the  system  becomes 

probabilistic. 

The  inputs  and  outputs  of  most  hydrologic  systems  do  not  interact 
at  a  single  point.   As  such,  there  should  properly  be  some  spatial 
description  of  the  system  processes.   The  complex,  dynamic,  nonlinear 
models  of  channel  routing  of  surface  runoff  typically  account  for  some 
of  the  spatial  variation  of  parameters  over  the  channel  or  catchment. 
The  parameters  of  simple  linear  models  usually  apply  over  the  entire 
channel  or  catchment  area  and,  therefore,  such  models  operate  as  if  all 
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input/output  interactions  occurred  at  a  single  point  (i.e.,  they  are 
lumped  models) . 

In  general,  a  system  can  have  several  inputs  and  outputs  as  shown 
in  Figure  2.1.   The  hydrologic  systems  analyzed  in  this  work  are  assumed 
to  have  a  single  input  and  a  single  output. 

A  causal  system  is  one  in  which  outputs  cannot  precede  inputs. 
All  physical  systems  are  causal.   This  restriction  limits  the  range  over 
which  integral  equations  describing  the  system  are  applied. 

The  stationarity  of  a  system  is  defined  by  whether  or  not  the 
input/output  relationship  changes  with  time.   If  a  system  is  stationary, 
identical  input  produces  identical  output. 

Linearity  is  the  property  that  opens  the  doors  to  the  storehouse 
of  mathematical  tools.   If  system  inputs  are  defined  as  xlf  system 
outputs  as  yt>  and  the  transformation  of  x  into  y  by  the  system  is 

symbolized  as  >  ,  then  given 

_>  y    and    x  >  y 


11  2        2 

a  system  is  said  to  be  linear  if 

a  .  x  +  a  •  x  >  a  •  y  +  a  •  y  (2>1) 

11     2   2         112   2 

where    a   and  a2  are  constants  (Thomas  1969,  p.  140). 

The  analysis  of  nonlinear  systems  is  much  more  difficult  than  that 

of  linear  systems  because  the  principle  of  superposition  applies  to 

linear  systems  (Thomas  1969,  p.  139).   Superposition,  which  follows 

directly  from  the  definition  of  linearity,  says  that  if  inputs  are  added 

or  scaled,  outputs  are  similarly  added  or  scaled.   Therefore,  the  system 

behavior  can  be  analyzed  independently  of  the  inputs.   This  is  a  point 

of  great  significance:   if  the  system  output  for  an  instantaneous  unit 
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input  can  be  found,  the  system  behavior  for  any_  input  can  be  obtained  by 
considering  the  actual  input  as  a  series  of  scaled  unit  inputs  at 
various  instants  of  time  (Liu  and  Liu  1975,  pp.  135-146).   Given  the 
system  output  for  an  instantaneous  unit  input  (called  the  impulse 
response),  the  output  for  a  stationary  linear  system  can  be  obtained  for 
any  input  with  the  convolution  integral 


,(t)  =    /  h(t-T)-x(T)  dT 


(2.2) 


where    x('),  y(* )  are  the  system  inputs  and  outputs  respectively,  and 
h(» )  is  the  system  impulse  response  (Johnson  and  Johnson  1975, 
pp.  137-140). 
Stationary  nonlinear  behavior  can  sometimes  be  simulated  by  a 
nonstationary  linear  model.   This  is,  in  fact,  exactly  the  strategy 
behind  the  nonlinear  state  space  model  developed  in  Chapter  IV.   The 
theoretical  basis  is  most  clearly  seen  by  recognizing  that  for  a  linear 
model,  the  input/output  relationship  can  be  represented  by  matrix 
multiplication  as 

'  I  ] 


1,1,1, 

12    3 


h   h   h 

11   12   13 

0   h   h 

22   23 

0   0   h 

33 


0   0    0 


1" 


(2.3) 


=  [o  ,  0  ,  0  ,  •  •  •  o J 

12    3  n 

where    Ii  =  the  input  at  time  i,  an  element  of  row  vector  U 

0.  =  the  output  at  time  j,  an  element  of  row  vector  0,  and 
h±   j  =  an  element  of  the  input/output  relationship  matrix _H. 
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Note  that  I  and  0  are  time  vectors,  in  contrast  to  the  state  vector  X.  of 

equation  1.1.   Matrix  _H  transforms  I_  into  0.   What  we  want  at  any  time, 

♦-..„   n   in  vector  0.   To  do  this  we  need  the 
t,  is  to  generate  one  entry,  0t»  LU   veuut  1L' 

appropriate  weights  to  apply  to  all  inputs  1^  Up  to  and  including  lt. 
These  weights  will  fill  one  column  of  matrix  H.   However,  it  is 
physically  impossible  to  observe  the  weighting  function  of  a  system.   An 
approach  which  circumvents  that  problem  is  to  generate  the  model  impulse 
response.   This  is  equivalent  to  solving  f or  0_  in  equation  2.3  when  I_  is 
all  zeroes  except  for  a  single  entry  of  1  in,  say,  the  1th  location. 
Call  this  vector  6_±.      The  output,  0_,  generated  will  be  the  ic  _row  of 
matrix  H.   The  input,  6^  selects  the  ltb  row  of  H  and  presents  it  as 
the  output,  0.   For  a  linear  nonstationary  system  the  rows  of  H_ will,  in 
general,  be  different  from  each  other  (i.e.,  each  row,  i,  will  represent 
the  system  impulse  response  at  time  i).   The  assumption  implicit  when 
using  a  single  system  impulse  response  in  the  convolution  integral  is 
that  the  system  is  linear  and  stationary.   In  this  case,  the  ith  row  of 
H  will  be  equal  to  the  (i-l)th  row  of  H_ shifted  one  position  to  the 
right.   Therefore,  equation  2.3  can  be  rewritten  for  a  stationary  linear 
system  as 


h 
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(2.4) 
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The  desired  weighting  functions  (i.e.,  the  columns  of  H)  are  seen  in 
equation  2.4  to  be  the  time-reversed  images  of  the  system  impulse 
response  for  a  stationary  linear  system.   Figure  2.2  shows  how  the 
convolution  integral  (equation  2.2)  exploits  this  to  compute  the  system 

output. 

Nonlinear  system  behavior  can  be  modelled  by  letting  the  columns 
of  H  in  equation  2.3  vary  each  time  step.   The  elements  on  any  diagonal 
will  not  be  equal,  as  is  the  case  for  the  linear  stationary  system 
represented  by  equation  2.4.   Therefore,  the  system  impulse  response 
(the  rows  of  H)  will  change  each  time  step,  thus  permitting  a 
nonstationary  linear  model  to  account  for  stationary  nonlinear 
behavior.   The  state  space  model  developed  in  Chapter  IV  emulates  the 
nonlinear  system  behavior  by  varying  the  weighting  function  based  on  the 
level  of  flow. 

A  General  Linear  Differential  Equation 

Linear  time-invariant  equations  can  be  used  to  approximate  the 
behavior  of  many  physical  systems.   This  is  a  good  approximation  when 
system  characteristics  change  very  slowly  relative  to  variations  of  the 
inputs  (Ogata  1967,  p.  307).   Channel  and  catchment  systems  can  be  well 
represented  by  linear  time-invariant  equations. 

A  general  differential  equation  for  a  nonstationary  linear  system 
can  be  written  as 


(t).^l2+a   l(t).il^l+  •  •  • 
,  n     n-1      j,-11-1 

dt  dt  (2.5) 


+  ai(t).l4r-+  aQ(t)-y(t)  =  x(t) 
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Figure  2.2   Calculation  of  Channel  Outflow  with 

the  System  Impulse  Response  Function 
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where    x(t)  =  the  system  input, 
y(t)  =  the  system  output, 

a0(t),  ai(t),  .  .  .,  an-i(t),  a^O  are  coefficients,  and 
the  initial  conditions  are  specified  for 

.n-1. 


y( 


O,  &&K    ....  and 


•y(t0) 


dt 


dt 


n-1 


The  system  is  linear  because  there  are  no  powers  or  products  of  y. 
Equation  2.5  can  be  rewritten  in  matrix  form  as 


djY(t), 
dt 


-  A(t)«Y(t)  =  B(t)-X(t) 


(2.6) 


where 


Y(t)  = 


y(t) 

dy(t) 
dt 


dn"Vt) 


dt 


n-1 


A(t)  = 


0 

0 
a0(t) 
"aTO 


a  <t) 


n-1 


(t) 


JtT 


B(t)  = 


JtT 


,  and 
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X(t)  = 


0 
x(t) 


If  the  system  is  stationary  the  coefficients  A_(t)  and  _B(t)  of  the 
underlying  differential  equation  2.5  are  not  functions  of  time.   In  this 
case  equation  2.6  can  be  rewritten  as 

fLlKt))  _  A.Y(t)  =  B.x(t)  (2.7) 

dt         

which  can  be  solved  using  the  integrating  factor  e  -  to  yield 
Y(t)  =e(t-t0)^Y(t0)  +  J  e(t-T)^B.X(x) 


dx 


(2.8) 


(Ogata  1967,  p.  315). 

The  first  term  on  the  right  hand  side  of  equation  2.8  is  the  zero 
input  response  of  the  system.  With  no  input  the  system  would  move 
toward  equilibrium  starting  from  this  point  when  t  =  tQ  .   The  second 
term  is  the  zero  state  response,  which  describes  the  motion  of  the 
system  driven  by  input  X(x  )  for  tQ  <^t    <_  t . 


Transformation  of  Moments  by  a  Stationary  Linear  System 

One  way  to  describe  a  system's  input  and  output  signals  is  to 
study  their  moments.   The  following  section,  on  the  use  of  transform 
methods  to  study  moments  of  stationary  linear  systems,  is  based  on  work 
by  Dooge  (1973,  pp.  113-115). 

The  Rth  moment  of  a  function  about  the  time  origin  is 
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u:(f)   =  /  f(t)-tR  dt  (2.9) 

and  the  Rth  moment  about  the  center  of  mass  Is 

Un(f)  =  /  f(t)-(t-U')R  dt  (2.10) 

R         -oo  1 

The  output  signal  of  a  stationary  linear  system  can  be  expressed 
with  the  convolution  integral  (equation  2.2)  in  terms  of  the  input 
signal  and  the  system  impulse  response.   The  relationships  between 
moments  of  the  input,  output  and  impulse  response  are  most  easily  seen 
by  transforming  equation  2.2  with  either  the  Fourier  or  Laplace 
transform.   Since  moments  about  the  center  of  mass  of  the  input  and 
output  signals  are  of  interest,  we  will  use  the  bilateral  Laplace 
transform  defined  by 

00 

FD(s)  =  /  f(t)-e"st  dt  (2.1D 

_oo 

where    f(t)  is  a  function  in  the  time  domain. 

Differentiating  equation  2.11  R  times  and  setting  s  =  0  produces 
the  following  expression  for  the  Rth  moment  about  the  origin  in  terms  of 
the  Laplace  transform  (Dooge  1973,  p.  114). 

„.<„  ,  h,1.A,w!h  (7-l2) 

ds 
The  convolution  integral  can  be  expressed  in  the  transform  domain  by 

taking  the  bilateral  Laplace  transform  of  equation  2.2  to  yield 

YB(s)  =  HB(s)-XB(s)  (2-13> 

where    YB(s),  HB(s)  and  XB(s)  are  the  Laplace  transforms  of  y(t),  h(t) 
and  x(t),  respectively  (Dooge  1973,  p.  114). 
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The  Rth  moment  of  the  output  about  the  origin  follows  from  equation  2.12 


„.(y)  ,  c-»R.-4[V>U  (2-14) 

ds 
Substituting  equation  2.13  for  YB(S)  gives 

U'Cy)  =  (-i)R-4tHB(8)-XB(8)ls-0  (2'15) 

ds 

Solving  equation  2.15  with  Leibnitz's  rule  for  continued  differentiation 

of  a  product,  and  substituting  for  the  moments  with  equation  2.12  yield 

R 


k=0 


where    fR)  indicates  a  combination  of  R  things  taken  k  at  a  time,  or 

*■  ky 

rRl  =       ,R! 

lkJ    k!-  (R-k)!   ' 

U'   (h)  =   the  (R-k)th  moment  about  the  origin  of  the  system 
R-k 

impulse  response,  and 

U'(x)  =  the  kth  moment  about  the  origin  of  the  system  input 
k 

signal  (Dooge  1973,  p.  114). 
If,  instead  of  about  the  origin,  the  moments  are  taken  about  the 
centers  of  mass  of  the  input,  output  and  system  impulse  response 
functions,  the  following  relationship  holds  (Dooge  1973,  p.  115) 

k=0 
Only  the  first  moment  about  the  origin  and  the  second  moment  about  the 
center  of  mass  are  needed  for  the  ensuing  work.   For  R  =  1  equation  2.16 
reduces  to 

U'(y)  =  U'(h)  +  U'(x)  (2.18) 

1       1       1 
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Since  U  (• )  =  0,  for  R  =  2  equation  2.17  becomes 


U  (y)  =  U  (h)  +  U  (x)  (2.19) 

2        2        2 

Therefore,  for  the  first  moment  about  the  origin  (the  time  lag)  and  the 

second  moment  about  the  center  of  mass  (the  variance),  the  moment  of  the 

system  output  is  the  sum  of  the  moments  of  the  system  impulse  response 

and  the  system  input  (Dooge  1973,  p.  115). 

Summary 

The  impulse  response  is  the  output  generated  by  an  instantaneous 
unit  input  to  a  system.   Because  the  impulse  response  function 
completely  describes  the  behavior  of  a  stationary  linear  system,  we 
would  like  to  determine  the  impulse  response  of  the  channel  or  catchment 
physical  systems.   The  best  description  of  the  real  world  that  we  have 
is  the  primary  model.   The  goal  of  linear  analysis  of  the  system  is  to 
determine  the  impulse  response  of  the  primary  model  at  various  reference 
flow  levels.   The  changing  properties  of  the  impulse  response  function 
with  flow  level  lend  insight  into  the  nonlinear  behavior  of  the  physical 

system. 

An  operational  way  to  determine  the  impulse  response  of  the 
primary  model  is  to  pulse  the  model  with  a  unit  input  and  observe  the 
output.   That  is,  in  general,  the  approach  that  will  be  taken  to  find  a 
system's  impulse  response.   An  analytical  solution  can  be  found  for  the 
impulse  response  of  the  Saint-Venant  equations  of  unsteady  flow  for 
downstream  waves  in  a  prismatic  channel.   To  derive  the  impulse 
response,  the  Saint-Venant  equations  must  be  linearized  about  a 
reference  discharge.   This  process  is  described  in  the  next  section. 
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Linearization  of  the  Salnt-Venant  Equations 
for  a  General  Prismatic  Channel 

To  determine  analytically  the  impulse  response  for  routing 
unsteady  flow,  the  Saint-Venant  equations  must  be  linearized  in  terms  of 
small  perturbations  about  a  reference  flow  condition.   The  steps 
presented  in  this  section  are  outlined  in  Figure  2.3.   Dooge  (1973,  pp. 
245-253)  presented  the  special  case  for  linearizing  the  Saint-Venant 
equations  for  a  rectangular  channel  and  Chezy  friction  formula.   The 
ensuing  derivation  follows  Schaake  (1980)  who  extended  Dooge's  work  for 
a  general  prismatic  channel  and  resistance  formula. 

The  Saint-Venant  equations  for  unsteady  flow  can  be  expressed  as 

d-A  +  ll=   B-ll+i^l-  q  (2-20) 

3t   <Tx    3T    31T^ 

and 

L3v   3y  +  v,8v  =_S^(U  -v)  (2.21) 

g  3t   3x   g3x    f    o    g'A  ^  x 

where    A  =  cross  sectional  area, 
Q  =  discharge, 
t  =  time, 

x  =  distance  along  the  channel, 
v  =  average  velocity, 
y  =  depth  of  flow, 
B  =  surface  width  of  flow, 
q  =  lateral  inflow  per  unit  width, 

u   =  the  x  component  of  the  velocity  of  the  lateral  inflow, 
g  =  the  acceleration  of  gravity, 
gf  =  the  friction  slope,  and 
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<;  =  the  channel  bottom  slope  (Henderson  1966,  p.  287). 
Equation  2.20  is  the  continuity  equation  for  a  channel.   The  effects  of 
the  forces  acting  on  a  control  volume  within  the  channel  are  modelled  by 
equation  2.21,  the  momentum  equation.   Since  the  continuity  equation  is 
already  linear,  only  the  momentum  equation  requires  linearization. 

Equation  2.21  is  first  rewritten  in  terms  of  the  dependent 
variables  Q  and  A.   To  this  end,  the  friction  term  and  channel  geometry 
are  expressed  as 

v  =  Q=k.(|)P./s7  (2.22) 

and 

B  =  B-Ar  (2.23) 

0 

where    B  is  the  top  width, 

BQ  is  the  top  width  for  flow  Qq »  and 

k,  p  and  r  are  described  below. 
Equation  2.22  is  a  general  relationship  for  friction  in  open  channel 
flow  which  can  be  related  to  the  well-known  Manning  or  Chezy  friction 
formula.   With  the  approximation 

R      A  •   A  (2.24) 

Rh     P     B 

where    R^  =  the  hydraulic  radius  and 

P  =  the  wetted  perimeter. 
Equation  2.22  becomes  Chezy' s  friction  formula  with 

k   =   C 

a  ,  (2.25) 

and 

p  =   1/2 

or  Manning's  friction  formula  with 


29 


k     =      1.49/n 
and  <2'26> 

p  =   2/3 
where    C  =  the  Chezy  flow  resistance  factor  and 

n  =  the  Manning  flow  resistence  factor  (Daily  and  Harleman 
1966,  pp.  297-298). 
Equation  2.23  can  approximate  most  regular  cross  sectional  shapes.   For 
a  rectangular  channel  r  =  0;  while  for  a  triangular  channel  r  =  1/2. 
However,  a  trapezoidal  channel  is  not  well  modelled  by  equation  2.23  for 
all  flow  ranges.   For  shallow  flow  a  trapezoidal  channel  behaves 
similarly  to  a  rectangular  channel;  while  for  deep  flow  it  can  be 
approximated  by  a  triangular  channel.   There  is  a  transition  flow  range 
where  equation  2.23  does  not  represent  the  change  in  top  width  with  area 
for  a  trapezoidal  channel. 
Substituting  2.23  into  2.22  and  rearranging  gives 


,   =       0 

5f   k2  .A2lp(l-r)+lJ 


(2.27) 


Let 

m  =  p(l-r)+l 
and 

Ws~ 


(2.28) 


(2.29) 


0 
where    m  is  analogous  to  the  dimensionless  kinematic  wave  parameter 

defined  by  Eagleson  (1970,  p.  250). 
Substitution  of  2.28  and  2.29  into  2.27  gives 
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(2.30) 
Substituting  equations  2.20  and  2.30  into  the  momentum  equation  gives 

(g.AW,.^+2.B.Q.A.£tB.A^  = 

(2.31) 

3+r      g'B  2   3+r-2m 
g.  b-  A   •  S   -  °-=>  Q  •  A 

0    ct2 
where  henceforth  the  terms  involving  q  will  be  ignored  because  we  are 
studying  prismatic  channels  with  no  lateral  inflows.   Finally, 
differentiating  2.20  with  respect  to  x,  differentiating  2.31  with 
respect  to  t  while  holding  the  left-hand  coefficients  constant,  and 

equating  the  common  ,,  \%      term,  results  in 


A    jn  32Q   „    32Q    32Q  _ 

—  -  v*  J  • -  2*  v  5 — 5—  -  = 

1  ax2       3x«3t   3t2 

(2.32) 

„  3Q   2'g*S0  3Q 
2«m«g'S  -^-  +  ^^t 

6  J  i  x      v   3 1 


Linearizing  equation  2.32  about  Q   A,  ,  vQ  and  BQ  forms  the  linear, 
second  order,  hyperbolic  equation 


g-A    o^  32Q   .     320    32Q 
'  B     0j  9x2       0  3x-3t    3t2 

°  (2.33) 

?•  ?•  <? 

2  m  g  S0 3x  +  ~ 3t 

0 


where    Aq ,  vQ  and  BQ  are  the  cross  sectional  area,  velocity  and  top 

width,  respectively,  for  the  reference  flow,  Q0  ( Schaake  1980), 
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The  impulse  response  of  the  linearized  form  of  equation  2.33, 
presented  in  Appendix  A,  is  a  complex  function  containing  a  Dirac  delta 
and  a  modified  Bessel  function  (Harley  1967,  p.  18).   A  simpler, 
approximate  solution  can  be  found  with  an  assumption  from  kinematic  wave 
theory.   The  3x«3t  and  3 t2  terms  of  equation  2.33  can  be  replaced  by  a 
Sx2  term  using  the  relationship 

1Q  _  _r.30-  (2.34) 

3t  "    3~x 


where 

dx 
c  =  -1—   (the  wave  celerity), 
dt 

This  substitution  is  routinely  made  in  the  derivation  of  diffusion 
routing  equations  (Koussis  1976).   Here  its  use  results  in  the 
simplified  linear  approximation  to  the  complete  equations 

9Q  ,  -  .3Q  -    Q0   .h^.BJO  (2.35) 

3t   C0  5¥   2.Bq.So        9x2 

where 

$2  =  F2- [l-m(2-m)l  (2,36) 

0 

in  which 

Q  /A 
F   =    °   °   (the  Froude  number  at  Q  )  (2.37) 

0    /g-A  /B  ° 

0   0 

and 

c   =  m«v   (the  kinematic  wave  velocity)  (Schaake  1980).      (2.38) 

0       0 

The  diffusion  form  of  equation  2.35  appears  often  in  channel  and 
pollutant  transport  modelling  (Fischer  1967,  Carslaw  and  Jaeger  1959, 
Bansal  1971).   Koussis  (1976)  followed  a  derivation  similar  to  the  one 
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presented  above,  but  was  interested  in  a  simpler  numerical  solution  to 

the  Saint-Venant  equations,  not  an  analytical  one.   The  interest  in  this 

section  is  in  an  analytical  solution  to  equation  2.35.   This  will 

provide  the  theoretical  basis  for  comparison  of  results  from  a  numerical 

solution  of  the  full  Saint-Venant  equations,  and  the  linear  and 

nonlinear  state  space  models  developed  later  in  this  chapter  and  in 

Chapter  IV,  respectively. 

The  influence  of  all  terms  in  the  Saint-Venant  equations  has  been 

retained  in  the  parabolic  form  of  equation  2.35.   The  impulse  response 

function  of  equation  2.35  is 

(x-c  t)2 
h(x,t)  =     X    >exP{-  -TT&-}  (2.39) 

4«ifD«  t3 

where 

D-  ^_-fl-S2)  (2'4°) 

2'B  «S   l     ; 

0   0 

(Harley  1966,  Schaake  1980).   The  impulse  response,  h(x,t),  defines  the 
flow  for  all  positive  x  and  t.   The  area  under  h(x,t)  is  unity, 

(i.e.,  /  h(x,t)  dt  =  1)  to  preserve  continuity,  and  the  units  are 

0 
time-1.   Equation  2.39,  the  solution  for  an  instantaneous  input  at 

x  =  0,  t  =  0  in  equation  2.33,  appears  in  many  scientific  disciplines  as 

a  solution  to  the  advection-dif fusion  equation  (Carslaw  and  Jaeger  1959, 

p.  50).   Since  the  diffusivity,  0,  is  constant  for  any  reference  flow, 

Q„ ,  Fickian  diffusion  is  represented  by  equation  2.39  (Slade  1968, 

pp.  80-81). 
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Insights  from  Linear  Systems  Theory 


Introduction 


The  analytically  derived  impulse  response  for  the  Saint-Venant 
equations  is  strictly  correct  only  for  small  variations  about  the 
reference  flow.   In  the  spirit  of  linear  approximations  for  nonlinear 
systems,  the  impulse  response  can  be  used  to  approximate  the  channel 
response  for  all  flow  ranges,  as  is  shown  in  Chapter  III. 

Equation  2.39  completely  describes  the  response  of  the  general 
prismatic  channel  for  flows  close  to  Q0 .   This  analytic  solution  can  be 
manipulated  to  gain  insights  into  the  behavior  of  the  Saint-Venant 
equations,  from  which  2.39  was  derived. 

The  Gaussian  form  of  the  solution  of  equation  2.39  for  fixed  t  is 
well  known  (Crank  1956,  pp.  9-11).   However,  Schaake  (1980)  recognized 
that  for  fixed  x  the  form  of  equation  2.39  is  an  inverse  Gaussian 
probability  density  function  (pdf )  .   Given  this  insight,  a  wealth  of 
information  can  be  obtained  from  previous  studies  of  this  statistical 
distribution  (Johnson  and  Kotz  1970,  pp.  137-153). 

The  Inverse  Gaussian  PDF  as  an  Analytical  Solution  of  the  Unsteady  Flow 
Equations 

Since  1871,  when  Barre  De  Saint-Venant  proposed  the  equations  of 
unsteady  flow  in  an  open  channel,  people  have  sought  an  analytical 
solution  for  equations  2.20  and  2.21.   No  analytical  solution  to  the 
full  Saint-Venant  equations  has  been  found.   An  analytical  expression 
for  the  outflow  hydrograph  would  be  ideal  for  studying  the  behavior  of  a 
channel.   However,  the  ideal  has  not  been  attained  because,  in  addition 
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to  the  lack  of  an  analytical  solution  of  the  Saint-Venant  equations, 
there  is  no  analytical  expression  for  a  general  inflow  hydrograph. 

There  is,  however,  an  analytical  solution  to  the  equations  of 
unsteady  flow  for  a  special  class  of  inflow  hydrographs  in  a  prismatic 
channel.   Because  equation  2.39,  the  equation  for  an  inverse  Gaussian 
pdf,  is  the  response  for  an  impulse  input  to  a  prismatic  channel,  an 
inverse  Gaussian  pdf  can  be  observed  at  any  point  along  the  channel. 
Parameters  of  the  inverse  Gaussian  pdf  are  determined  by  the  channel 
properties  and  the  distance,  x,  from  the  point  at  which  the  impulse  is 
input.   As  shown  in  Figure  2.4,  an  inverse  Gaussian  pdf  with  x  =  xQ +  l 
is  the  analytical  expression  for  the  outflow  hydrograph  when  an  inflow 
hydrograph  is  specified  by  an  inverse  Gaussian  pdf  with  x  =  xQ .   Channel 
properties  determine  the  parameters  Cq  and  D  of  equation  2.39.   This 
means  that  for  the  class  of  inflow  hydrographs  defined  by  inverse 
Gaussian  pdfs,  the  outflow  hydrographs  from  a  prismatic  channel  are  also 
inverse  Gaussian  pdfs  (Schaake  1980)!   This  is  of  great  significance 
because  the  parameters  of  the  outflow  inverse  Gaussian  pdf  are  related 
simply  to  the  parameters  of  the  inflow  inverse  Gaussian  pdf  and  the 
channel  properties. 

Moments  of  the  Inverse  Gaussian  PDF 

The  time  delay  and  dispersion  characteristics  of  2.39  are 
expressed  by  its  first  two  moments  (Johnson  and  Kotz  1970,  pp.  137- 
140).   The  first  moment  about  the  origin  (i.e.,  the  mean)  describes  the 
time  lag  properties  as 


k   =  t.  =  

i     I        c 
0 


(2.41) 
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The  second  moment  about  the  mean  (i.e.,  the  variance)  describes  the 
dispersion  properties  as 


Q  «x 

k  =  a2  =  S >(l-*' 

2         B  •  S  •  C3 

0   0   0 


(2.42) 


The  coefficient  of  variation,  a  dimensionless  measure  of  the  impulse 
response,  can  be  obtained  from  2.41  and  2.42  as 

1/2 


0 
0 
B  • S  •  c  «x 

0   0   0 


1-4' 


(2.43) 


The  form  of  the  coefficient  of  variation,  equation  2.43,  indicates 
that  the  kinematic  wave  assumption  is  most  nearly  correct  for  shallow 


U 
flow  (small   <°-  )  on  steep  slopes  (large  SQ ) .  Woolhiser  and  Liggett 

0   0 


■Q  ,.     ..uunu,,,,   .x,,.„   ^gget 

(1967)  define  a  dimensionless  parameter  to  determine  when  the  kinematic 


wave  assumption  is  appropriate 


(2.44) 


L  -S 

k   =  _0 0_ 

WL   H  .p2 

0   0 


where    L   =  characteristic  length, 
Sn  =  channel  bottom  slope, 
Hq  =  depth  of  flow,  and 

F2  =  Froude  number. 

0 

If  kWL  >  10,  the  Saint-Venant  equations  are  approximated  well  by  a 
kinematic  wave  solution  (Eagleson  1970,  p.  336).   These  conditions  are 
generally  met  by  overland  flow,  with  its  shallow  depth  of  flow  and  steep 
slopes.   For  river  channels  with  steep  slopes,  the  SQ  term  of  the 
momentum  equation  dominates,  making  the  kinematic  wave  assumption  a  good 
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approximation  of  the  Saint-Venant  equations.   For  rivers  with  flat 
slopes,  the  S  and  ^-   terms  are  dominant  in  equation  2.21  (Henderson 
1966,  p.  364).   The  diffusion  analogy  model,  discussed  in  Appendix  A, 
uses  only  these  two  terms  for  channel  routing. 

The  Inverse  Gaussian  PDF  as  an  Inflow  Hydrograph 

The  interrelationships  between  peak  flow  (Q_),  celerity  (cQ)  and 
timing  properties  (as  measured  by  the  variance,  o2^)  of  the  inflow 
hydrograph  can  be  expressed  for  a  given  total  flow  volume  (V) 
(Schaake  1980).  With  these  relationships  several  properties  of 
hydrographs  input  to  routing  models  can  be  controlled.   Any  three  of  Qpi 

c   a2   or  V  can  be  specified  arbitrarilv.   Since  equation  2.39  is  the 

0  '   in 

solution  for  a  pulse  input  at  x  =  0  and  t  =  0,  we  will  define  a  point 
downstream  at  x  =  Xq  to  observe  the  flow  for  t  >_  0.  With  this  notation 
and  h(«,-)  defined  by  equation  2.39,  the  total  flow  volume,  Inflow 
hydrograph,  and  peak  discharge  are,  respectively, 


V  =  f  Q(x  ,t)  dt  =  /  V«h(x  ,t)  dt  =  V-J  h(x  ,t)  dt 
0  0  o 


(2.45) 


Q(x  ,t)  =  Vh(x  ,t)  (2-46) 

0  o 


Q  =  Vhfx  ,t  (x  ), 

yp      v  0   P   0 


where  h(x,t)  is  a  maximum  when 


(2.47) 


tp(,>-i£i{i*^C)2i'"  -  >  i 


3«D; 
c   '  — 
0 
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By  choosing  any  three  of  V,  Qp,  a\n   and  Cq ,  a  hydrograph  with 
certain  desired  properties  can  be  selected.   Since  equation  2.39  defines 
the  response  to  an  impulse  input,  an  inverse  Gaussian  shape  hydrograph 
can  be  observed  at  any  point  along  the  channel.   The  channel  can  be 
divided  into  two  sections  as  shown  in  Figure  2.4.   The  upstream, 
ficticious  section  has  an  impulse  input  and  produces  an  inverse  Gaussian 
pdf  outflow  hydrograph  which  becomes  the  inflow  hydrograph  for  the 
downstream,  real  channel  of  interest.  When  an  impulse  is  input  to  the 
upstream  section  at  x  =  0,  an  inverse  Gaussian  pdf  appears  as  input  to 
the  downstream  section  at  x  =  Xq .   The  variance  of  the  flow  transient  as 
x  =  xQ  (a2,  )  can  be  selected  by  proper  choice  of  Xq  .   Letting  a2  and  x 
of  equation  2.42  equal  a2      and  xQ ,  respectively,  gives 

a2   -   Q0  X0  -h^2]  (2.49) 

0   0   0 
Substituting  equation  2.40  into  equation  2.49  simplifies  the  expression 
for  the  inflow  variance  to 


(2.50) 


The  parameters  D  and  c.  are  specified  by  the  channel  geometry  and 
slope.   The  value  of  xq  can  be  arbitrarily  selected  to  give  the  desired 
inflow  variance.   Equation  2.50  can  be  rearranged  to  solve  for  xQ  as 


c3-o2 

,  =  _0 In  (2.51) 

0      2-D 


2 

2'D'x 

0 

in 

c2 
0 
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The  channel  section  where  0  <_  x  <_  Kq  is  a  ficticious  reach  designed 
solely  to  generate  an  inflow  hydrograph  with  desired  properties  for  the 
real  channel  being  modelled.   This  real  channel  is  defined  by  the 
x  _<_  x  <_xQ  +  L  channel  section  where  L  is  the  length  of  the  real 
channel.   The  outflow  hydrograph  is  also  an  inverse  Gaussian  shape  with 
variance 

2'D'fx  +  L] 


a' 


_0 (2.52) 


out        c3 

0 

Equations  2.18  and  2.19  can  be  rewritten  to  solve  for  the  mean  and 
variance  of  the  impulse  response 

U'(h)  =  U'(y)  -  U'(x)  (2.53) 

1  1       1 

and 

U  (h)  =  U  (y)  -  U  (x)  (2.54) 

2  2        2 

With  these  relationships  and  equations  2.41  and  2.50,  the  mean  of  the 

impulse  response  function  for  the  linearized  Saint-Venant  equations  can 

be  expressed  as 

x  +  L   x 
u»(h)  =_Q JL  (2.55) 

1       c    c 

1  0      0 


U'(h)  =  lagK    .  -i-  (2.56) 

i         channel   c 

1  0 


an 


d  the  variance  of  the  impulse  response  function  is 

2«D'(x  +  L)    2-D'x 

0(h)  -o2    -o]      - B. 2.  (2.57) 

2        out     in        c3  c3 

0  0 
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„  (h)  .0*     .    -*£*.  (2.58) 

o       channel      3 

0 
Equations  2.56  and  2.58  express  the  lag  and  variance,  respectively, 

which  are  added  to  an  input  hydrograph  as  it  passes  through  a  prismatic 

channel  specified  by  parameters  D,  cQ  and  L.   The  channel  impulse 

response  will  vary  depending  on  the  reference  flow  level  chosen,  since  D 

and  c  are  functions  of  Q0 .   Note,  however,  that  the  channel  impulse 

response  is  independent  of  xQ  .   This  means  that  once  Qq  is  chosen,  the 

impulse  response  of  the  system  is  not  dependent  on  the  input  signal. 

The  development  of  a  linear  state  space  model  which  follows  later 

in  this  chapter  uses  the  properties  of  the  linearized  Saint-Venant 

equations  described  above.   The  remainder  of  this  section  follows 

Schaake  (1980)  and  focuses  on  several  properties  of  equation  2.39  that 

lend  additional  insight  to  routing  in  hydrology. 

Additional  Properties  of  the  Inverse  Gaussian  PDF 

Attenuation  of  the  peak  inflow,  Q   to  the  real  section  of  the 
channel  can  also  be  determined  from  an  analysis  of  equation  2.39. 
Define  the  attenuation  of  the  peak  inflow  at  a  point  distance  L  below 
the  upstream  location,  Xq  »  as 

V'hfx  +  L,t  (x  +  L)' 


A(L)  = 


V'hfx  ,t  (x  ), 
0   P   0 


g  0  (2.59) 


where  the  numerator  is  the  peak  outflow  rate  at  location  xQ +  L  and  the 
denominator  is  0p-   Cancelling  the  common  factor  ¥  produces 

hfx  +  L,t  (x  +  L)) 
A(L)  =  3 L_J (2.60) 

h(x  ,  t  (x  )) 
0    P  0 
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Schaake  (1980)  notes  the  interesting  result  that  although  A,  the  inflow 
hydrograph  attenuation,  is  independent  of  V,  the  downstream  peak  outflow 
rate  is  directly  proportional  to  V. 

The  channel  lag  of  the  center  of  mass  (i.e.,  the  mean)  of  the 
inflow  hydrograph  is  given  by  equation  2.56.   The  peak  lag  time  can  be 
derived  from  equation  2.39  as 


C  „  =  t  fx  +  LJ  -  t  (x 
P*     P^  0         P   0 


1/2 


1/2 


'pi 


1  + 


c  «(x  +  L 

0    0 
3^D 


!  +  J  0   0 


"snr 


(2.61) 


(2.62) 


A  Linear  State  Space  Approach 


Introduction 


The  inverse  Gaussian  pdf  is  a  linear  approximation  of  the  full 
nonlinear  Satnt-Venant  equations.   For  small  variations  in  flow  the 
Inverse  Gaussian  pdf  produces  results  nearly  Identical  with  the  full 
equations  (Chapter  III).   Unfortunately,  derivation  of  a  state  space 
model  for  equation  2.39  is  cumbersome.   A  state  space  model  with 
parameters  related  to  the  parameters  of  the  inverse  Gaussian  pdf  Is 
derived  in  this  section. 

The  inverse  Gaussian  pdf  is  defined  by  equation  2.39  in  terras  of 
the  parameters  cQ ,  the  wave  velocitv,  and  0,  a  diffusion  constant 
(Johnson  and  Kotz  1970,  p.  137).   This  form  of  the  inverse  Gaussian  pdf 
follows  naturally  from  the  physically  based  parameters  of  the 
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Saint-Venant  equations.   Equations  2.38  and  2.40  define  Cq  and  D, 
respectively,  in  terms  of  the  reference  discharge  and  cross  sectional 
properties.   However,  in  the  ensuing  analysis,  the  system  impulse 
response  will  be  determined  from  its  moment  properties  which  will  be 
deduced  from  the  changing  moments  of  the  inflow  and  outflow,  using 
equations  2.56  and  2.58.   To  facilitate  the  determination  of  the  inverse 
Gaussian  pdf  from  its  moment  properties,  an  alternative  form  of  equation 
2.39  can  be  defined  in  terms  of  two  different  parameters:  u,  the  time 
mean  of  the  distribution,  and  X,  an  inverse  measure  of  dispersion 
(Johnson  and  Kotz  1970,  p.  138).   Equation  2.39  can  be  rewritten  in 
terms  of  \i   and  X  as 

h(x,t)H  fIG(t|u,X)  =/^^7^-exp{-X.(t-y)2/(2.y2.t)}    (2.63) 
The  parameters  of  equations  2.39  and  2.63  are  related  by 

u  =  JL  (2.64) 

c 

0 

k    _  x2  (2.65) 

The  mean  (y  IG)  f    variance  (Ojq),  and  skewness  (Y  iq)  of  the  inverse 
Gaussian  pdf  are  simple  functions  of  the  parameters  u  and  X  (Johnson  and 
Kotz  1970,  pp.  139-140). 

^IG 

02   =p3/X  (2.67) 

IG 

IG 
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These  simple  relationships  facilitate  the  determination  of  state  space 
model  parameters.   This  is  demonstrated  in  the  following  sections. 

Criteria  for  Matching  Impulse  Response  Functions 

The  value  of  a  linear  system  analysis  of  the  unsteady  flow 
equations  is  that  is  allows  the  development  of  a  simple  model  that  has 
an  impulse  response  function  similar  to  the  impulse  response  of  the  full 
Saint-Venant  equations.   For  small  perturbations  about  a  reference  flow, 
q  ,  equation  2.63  defines  a  system  impulse  response  which  is  the  same  as 
the  impulse  response  for  the  full  Saint-Venant  equations.   In  Chapter 
III,  we  will  see  that  equation  2.63  is  an  excellent  approximation  to  the 
Saint-Venant  equations  for  small  (10  percent  of  QQ  )  inflow  transients. 
The  linear  approximation  deteriorates  for  large  (2000  percent  of  Q0 ) 
inflow  hydrographs.   The  results  are  also  presented  in  Chapter  III  to 
help  show  the  limitations  of  the  linear  theory.   The  strategy  employed 
to  determine  the  limitations  of  the  linear  theory  is  depicted  in  Figure 
2.5.   The  primary  thesis  of  this  work  is  that  the  fundamental 
nonlinearities  present  when  routing  large  inflow  hydrographs  with  the 
Saint-Venant  equations  can  be  accounted  for  in  a  model  that  matches  the 
impulse  response  of  the  Saint-Venant  equations  at  all  reference  flow 

levels,  Q0 • 

There  are  a  number  of  possible  criteria  which  could  be  used  to 
compare  impulse  response  functions  for  the  simple  and  complete  models. 
In  general,  functional  relationships  for  the  impulse  responses  cannot  be 
obtained,  so  the  criteria  for  matching  impulse  response  functions  must 
be  able  to  compare  the  responses  of  the  systems  to  impulse  inputs.   Some 
of  the  methods  available  for  comparison  of  impulse  responses  are 
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1.  least  squares  regression, 

2.  maximum  likelihood  analysis, 

3.  spectral  analysis,  and 

4.  matching  moments. 

Direct  derivation  of  parameters  for  the  least  squares  criterion 
may  not  be  easy.   For  the  simplest  of  models,  a  functional  relationship 
may  be  found  between  the  parameters  and  the  criterion  to  be  minimized. 
Differentiation  of  that  function  for  each  parameter  and  solution  of  the 
set  of  the  simultaneous  equations  generated  is  almost  always  difficult 
because  of  nonlinearities  (Dooge  1970,  pp.  170-171).   Statistical 
packages  are  available  to  perform  least  squares  regression  and  to 
determine  parameter  values  but,  for  any  except  the  simplest  models, 
those  packages  can  be  costly  to  use. 

Maximum  likelihood  parameter  identification  suffers  from  the  same 
weaknesses.   It  may  be  difficult  to  determine  the  likelihood  function 
for  nonlinear  models  and  the  cost  of  using  statistical  packages  to 
identify  parameters  involves  significant  cost  for  complex  models. 
Spectral  analysis,  which  involves  transformations  into  the 
frequency  domain,  comparison  of  selected  properties,  and  transformation 
back  into  the  time  domain,  in  principle,  can  be  carried  out  to  determine 
parameter  values.   The  state  space  and  full  systems  can  be  equated  by 
matching  pole  and  zero  locations  of  the  transfer  functions  in  the 
transform  domain.   Since  most  hydrologic  situations  are  not  described 
functionally,  the  transformations  to  and  from  the  spectral  domain  must 
be  carried  out  numerically.   Numerical  inversion  of  the  transform  back 
to  the  time  domain  is  almost  always  difficult  (Dooge  1973,  p.  31).   The 
autocovariance  properties  of  the  impulse  response  function  previously 
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have  been  used  to  find  the  parameters  of  a  state  space  model  for  a  unit 
hydrograph  (Goldstein  and  Larimore  1980).   This  approach  was  not 
pursued,  but  may  prove  useful  for  determining  state  space  model 
parameters  when  the  impulse  response  is  not  a  single-peaked  function. 
The  mean  and  variance  of  the  impulse  response  functions  for  the 
linear  and  full  systems  can  be  compared  to  determine  parameter  values  by 
the  method  of  moments.   If  the  first  n  moments  of  two  impulse  responses 
are  identical,  the  two  systems  will  produce  identical  output  for 
polynomial  input  of  degree  n  or  less  (Dooge  1970,  pp.  170-171).   The 
method  of  moments  is  used  for  subsequent  analysis  of  the  state  space 
model  because  it  is  simple  and  allows  the  convenient  parameterization  of 
the  model  in  terms  of  the  moments  of  the  impulse  response  function. 

Computing  the  Parameters  of  a  Gamma  PDF 

A  cascade  of  linear  storage  reservoirs  is  a  well-studied  routing 
model  (Zoch  1934,  Nash  1959,  Dooge  1973,  p.  176).   This  model  is 
amenable  to  state  space  formulation  (Szollosi-Nagy  1981)  and  can  be 
functionally  represented  by  a  3-parameter  gamma  pdf  (Johnson  and  Kotz 
1970,  p.  166) 

k.((t-a)»k)n~1-exp{-(t-a)'k}  (2>69) 

c(t) ^ 

where    a  =  the  pure  delay, 

n  =  the  number  of  reservoirs, 

k  =  the  time  delay  in  each  reservoir,  and 

r(«)  =  the  gamma  function  (Abramowitz  and  Stegun  1964,  p.  255). 
As  formulated  in  equation  2.69,  n  can  take  on  any  positive 
value.   The  interpretation  of  n  as  the  number  of  reservoirs  restricts  it 
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to  positive  integers,  implying  that  T(n)  could  be  replaced  by  (n-1)! 
(the  factorial  function)  in  2.69.   Since  it  is  useful  to  work  with 
continuous  parameters,  the  positive  integer  restriction  on  n  will  be 
relaxed  temporarily.   The  concept  of  integer  values  for  n  will  be 
reinstated  for  implementation  of  the  state  space  model. 

Equations  for  the  mean  (MG),  variance  (Oq),  and  skewness  (y q)   of 
the  3-parameter  gamma  pdf  are 


UG 


a  +  n/k  (2.70) 

2  =  n/k2  (2.71) 


Y   =  2//k 
G 


(2.72) 


(Johnson  and  Kotz  1970,  pp.  167-168).   Parameters  of  the  inverse 
Gaussian  and  gamma  pdfs  can  be  related  by  equating  their  means, 
variances  and  skewnesses  as  defined  by  equations  2.66  through  2.68  and 
2.70  through  2.72.   When  the  gamma  pdf  parameters  are  found  in  terms  of 
the  inverse  Gaussian  parameters,  the  surprisingly  simple  results  are 

a-^.y  (2'73) 

=  4_   X_  (2.74) 

n   9   H 

k  -  2  .  L  (2-75) 

T   M2 

As  demonstrated  in  Figures  2.6  and  2.7,  an  inverse  Gaussian  pdf 
can  be  approximated  well  with  a  gamma  pdf  using  equations  2.73  through 
2.75.   The  relationships  based  on  matching  moments  begin  to  deteriorate 
as  n  approaches  1,  because  the  gamma  pdf  becomes  simply  a  shifted 
exponential  (Johnson  and  Kotz  1970,  p.  166).   Results  for  this  case  are 
shown  in  Figure  2.8. 
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Figure  2.6   Comparison  of  the  Inverse  Gaussian 
and  Gamma  PDFs  with  Parameters 
Matched  by  Method  of  Moments 
for  a  =  1,  n  =  4  and  k  =  2.0 
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Figure  2.7   Comparison  of  the  Inverse  Gaussian 
and  Gamma  PDFs  with  Parameters 
Matched  by  Method  of  Moments 
for  a  =  8,  n  =  8  and  k  =  0.5 
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Figure  2.8   Comparison  of  the  Inverse  Gaussian 
and  Gamma  PDFs  with  Parameters 
Matched  by  Method  of  Moments 
for  a  =  1,  n  =  1  and  k  =  0.5 
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Parameter  values  for  a  cascade  of  linear  reservoirs  model  can  be 
computed  from  values  of  y  and  X  with  equations  2.73  through  2.75. 
Values  of  M   and  X  can  be  determined  from  channel  properties  with 
equations  2.64  and  2.65  if  x,  c   and  D  are  known,  or  from  the  mean  and 
variance  of  a  channel  response  function  with 

(2.76) 


U3 
X    =  -11  (2.77) 

a2 
IG 


(Johnson  and  Kotz  1970,  pp.  139-140).   Equations  2.76  and  2.77  can  be 
substituted  into  equations  2.73  through  2.75  to  provide  relationships 
for  the  parameters  of  the  gamma  pdf  in  terms  of  the  mean  and  variance  of 
the  system  impulse  response  function  (i.e.,  the  inverse  Gaussian  pdf)  as 

a   T   WIG 

M2 
n  =  1.  .11  (2.79) 

v   ~  2  ■   IG  (2.80) 

3   a2 
IG 

In  practice  the  mean  (u_J  and  variance  (a2  )    would  be  estimated  from 
XG  •*■  J 

the  moments  of  the  primary  model  input  and  output  hydrographs .   To  be 
strictly  correct  in  this  case,  the  inverse  Gaussian  mean  and  variance  in 
equations  2.78  through  2.80  should  be  replaced  by  uZ1Q   and  a2^, 
respectively,  where  the  *  indicates  a  value  estimated  from  data. 

The  goal  of  this  approach  is  to  emulate  a  calibrated  primary 
routing  model  (such  as  the  numerical  solution  of  the  Saint-Venant 
equations)  with  a  linear  state  space  model.   The  general  approach  is  to 
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1.  simulate  with  the  primary  model  using  an  inverse  Gaussian 
pdf  as  the  input  hydrograph, 

2.  compute  the  mean  and  variance  of  the  input  and  output 
hydrographs, 

3.  estimate  the  mean  (ji-J  and  variance  [o\G)      of  the  system 
impulse  response  with  equation  2.56  and  2.58,  and 

4.  compute  a,  n  and  k  from  equations  2.78  through  2.80. 
When  the  values  of  a,  n  and  k  computed  in  step  4  above  are 

substituted  into  equation  2.69,  G(t)  defines  the  impulse  response  of  the 
desired  linear  system.   A  state  space  model  with  an  impulse  response 
identical  to  equation  2.69  is  derived  in  the  next  section. 

A  Cascade  of  Linear  Reservoirs  State  Space  Model 

The  state  space  model  developed  in  this  section  is  based 
conceptually  on  a  cascade  of  linear  reservoirs.   Figure  2.9  provides  an 
overview  of  the  steps  taken  to  develop  the  linear  state  space  model. 
The  relationship  among  inflow,  outflow  and  storage  for  a  reservoir  is 
described  by  the  lumped  continuity  equation 

g!-i-o  (2-81) 

dt 

where    S  =  the  reservoir  storage, 

I  =  the  reservoir  inflow,  and 

0  =  the  reservoir  outflow. 
A  reservoir  is  linear  when  the  relationship  between  storage  and  outflow 
is  linear  as 

S=<-0  <2'82) 

where    <  =  a  constant. 


53 


RESERVOIR 

CONTINUITY 

EQUATION 

AS/At  -  I  -  0 


LINEAR 

RESERVOIR 

EQUATION 

S  -  K-0 


GENERAL  STATE  SPACE 
EQUATION 


Xfi  -  E>Xt 


fclLai 


CHANNEL  IMPULSE 

RESPONSE  FUNCTION 

AT  REFERENCE 

FLOW  LEVEL 


EQUATIONS  FOR 

LINEAR  RESERVOIRS 

IN  A  CASCADE 


± 


EXPRESS  F  AND  G 
IN  TERMS  OF 
a,  n  AND  k 


COMPUTE  a,n,k 

AT  REFERENCE 

FLOW  LEVEL 


7 


STATE  SPACE  MODEL  WITH 

F  AND  G  AS  FUNCTIONS 

OF  CHANNEL"  IMPULSE  RESPONSE 


Figure  2.9    Development  of  the  Linear 
State  Space  Model 
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The  state  space  model  will  operate  with  discrete  time  steps  of  At 
when  implemented,  so  rewriting  equation  2.81  for  discrete  time  yields 

Ji-I-0  (2.83) 

At 

where    T  =  (1-P  )•  I   +  P'  I 
1      2 

0  =  (l-p)'O  +  P'O 

1       2 

Subscripts  1  and  2  indicate  the  beginning  and  end  of  a  time  step, 
respectively,  and  p  can  take  on  values  from  0  to  1 ;  I  and  0  are 
weighted  average  values  of  the  inflows  and  outflows,  respectively,  at 
the  beginning  and  end  of  a  time  step. 
Substituting  for  T  and  0  in  equation  2.83  gives 

7^-  (l-p)-I   +P-I   -  (l-p)-O  -p'O  (2.84) 

At  12  12 

Rewriting  AS  in  terms  of  storage  at  the  beginning  and  end  of  a  time  step 
and  substituting  for  S  with  equation  2.82  produce 

k«0  -  <*0 

2-_ L=(l-p).i  +p-I  -  (l-p).O  -P'O  (2.85) 

At  12  12 


K'O   +  P'O  'At  =  Af((l-p)'l   +  P'T  )  + 
2       2  1       2 


(2.86) 


K'O   -  (l-p)'O  «At 

1  1 


or,  solving  for  02 

(k  -  (l-p)-At) 

°2  -jrrhnr^1-**1!*''1!  +    mp^o  >0i        (2-87) 

Equation  2.87  applies  for  any  reservoir  in  a  cascade.   The  form  of 
equation  2.87  varies  slightly  for  the  first  reservoir  because  the 
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inflow,  It.  is  externally  specified  for  the  first  reservoir.   For  all 
subsequent  reservoirs,  the  inflow  is  equal  to  the  outflow  from  the 
previous  reservoir,  as 

I    =0       ;   i  >  1  (2-88) 

t,i    t,i-l 

where  the  first  subscript  indicates  the  time  step  and  the  second 
subscript  indicates  the  relative  position  of  the  reservoir  in  the 
cascade. 

Rewriting  equation  2.87  with  this  double  subscript  notation  for 

the  first  reservoir  yields 

[k    -   (l-p)-At) 

Vi.i  ■  T-ta-^t  +  »'hJ  + 1 '°t,i  (2'89) 

where   It  ±s   the  inflow  to  the  first  reservoir  at  time  t  and 

E,    =  k  +  p*At. 
The  governing  equation  for  all  subsequent  reservoirs  is 

Vi,i=^(H),0t,i-i+p,Vi,iJ  + 

(2.90) 


[k    -  (l-o)-At) 

1 '"t.i 

where    i  >  1 . 

Recall  the  general  form  of  a  state  space  model  from  equation  1.1 

X    =  F    «X   +  G    *U   , 
^t+1    -X,t  -t   -X,t  -t+1 

Equations  2.89  and  2.90  must  be  organized  into  the  matrix  form  indicated 
by  equation  1.1.   For  a  linear  system  the  state  transition  and  input 
coefficient  matrices  are  constant,  so  equation  1.1  can  be  simplified  to 

X    =  F'X  +  G'U   ,  (2.91) 

^t+1 1 1+1 
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Defining  the  state  vector  for  a  cascade  of  n  linear  reservoirs  as 


X.  = 


t,l 


t,2 


t,n 


(2.92) 


and  the  system  input  vector  as 


U-P)-it  +  P-it+1 


Setl 


(2.93) 


F  and  G_  can  be  determined  by  substituting  equations  2.89,  2.90,  2.92  and 

2.93  into  equation  2.91. 

Defining 

'<    -  (l-p).At) 


*  = 


T 


T  = 


At 


T 
a   =  (i-p)  +  p«* 


(2.94) 
(2.95) 

(2.96) 


the  F  and  G  matrices  for  a  cascade  of  n  identical  linear  reservoirs  are 


F  = 


T»n 


.2 


P'T^'ft 


n-2   n-1 
p    *T    'U 


0       0 

y  o 

t  •  n    *     o 


p«T  '71        T 


Q    * 


(2.97) 


and 
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G  = 


P'T' 


pn_1.Tn 


(2.98) 


Parameters  k  (equation  2.80)  and  <    (equation  2.94)  are  related  by 

k  =  1/k  (2.99) 

The  parameter  a  (equation  2.69)  is  accounted  for  in  the  state  space 
model  by  delaying  the  time  at  which  the  outflow  from  the  ntn  reservoir 
appears  by  a  units  of  time. 


Summary 

A  linear  state  space  model  that  emulates  the  behavior  of  the  full 
Saint-Venant  equations  about  a  reference  flow  level  has  been 
developed.   Parameters  of  the  state  space  model  are  determined  by 
matching  the  mean  and  variance  of  the  system  impulse  responses  of  the 
primary  and  surrogate  models.   Verification  of  the  cascade  of  linear 
reservoirs  state  space  model  by  comparing  the  results  with  equation  2.39 
(the  impulse  response  of  the  linearized  Saint-Venant  equations)  and  a 
study  of  the  limitations  of  the  linear  state  space  model  for  large 
inflow  hydrographs  are  the  major  topics  of  Chapter  III.   The  work  of 
Chapter  IV,  in  which  the  important  nonlinearities  are  analyzed  and  a 
nonlinear  state  space  model  is  developed,  is  motivated  by  the  -problems 
arising  when  large  inflow  hydrographs  are  routed  with  the  linear  state 
space  model. 


CHAPTER  III 
LIMITATIONS  OF  THE  LINEAR  THEORY 


Introduction 


The  linear  theory  developed  in  Chapter  II  for  prismatic  channels 
is  applicable  strictly  only  for  small  perturbations  around  a  reference 
flow  level.   In  this  chapter,  outflows  predicted  by  the  linearized 
Saint-Venant  equations  and  the  linear  state  space  model  are  compared 
with  numerical  solutions  of  the  full  Saint-Venant  equations  for  both 
small  (10  percent  of  baseflow)  and  large  (2000  percent  of  baseflow) 
inflow  hydrographs. 

The  model  used  as  the  numerical  solution  of  the  Saint-Venant 
equations  is  the  National  Weather  Service  Dynamic  Wave  Operational 
(DWOPER)  Model  (Fread  1978).   This  dynamic  wave  routing  model  is  based 
on  an  implicit  finite  difference  solution  of  the  complete 
one-dimensional  Saint-Venant  equations  of  unsteady  flow. 


8(A+A  ) 
9Q+_0   _q=0  (3.1) 

dX        d  t 


jt+-^__+g.A.(|i+Sf)-,.v^«f.B-0  (3.2, 

where  0  =  discharge, 

A  =  cross  sectional  area, 

A  =  off-channel  storage  area, 

q  =  lateral  inflow  or  outflow, 
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x  =  distance  along  the  channel, 

t  =  time, 

g  =  the  acceleration  of  gravity, 

h  =  the  water  surface  elevation, 

vx  =  velocity  of  lateral  inflow  in  the  x-direction, 

Wf  =  the  wind  term, 

B  =  channel  top  width,  and 

gf  =  the  friction  slope  defined  as 

«  -   "MOM  (3.3) 

sf 2 — S7T 

2.2' A  'R 
in  which   n  =  the  Manning  roughness  coefficient, 

R  =  the  hydraulic  radius,  and 

|«  |  =  the  absolute  value  function  (Fread  1978). 
Off-channel  storage,  lateral  inflows  and  wind  effects  are  ignored  in 
this  work. 

The  solution  technique  in  the  DWOPER  model  is  implicit;  therefore, 
the  time  step  size  can  be  selected  based  on  accuracy  rather  than 
numerical  stability.   This  makes  the  DWOPER  model  very  efficient  in  the 
use  of  computer  time.   Simulation  of  actual  river  reaches  typically 
requires  from  10  to  30  seconds  of  CPU  time  on  an  IBM  360/195  computer 
(Fread  1978). 

Verification  of  the  Linear  Saint-Venant  Equations 

The  impulse  response  function  of  the  linearized  Saint-Venant 

equations  was  derived  in  Chapter  II  as  equation  2.39. 

(x-c  t)2 
h(x,t)  ^expj-^ng-p — J 

WD' t3 
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where 


and 


Q 

2'B  «S   ^ 

0   0 


;   =  m*  v 

0       0 


In  this  section,  results  obtained  by  simulating  with  the  DWOPER  model 
and  by  solving  equation  2.39  are  presented  for  10  percent  of  baseflow 
transients  on  baseflows  of  10000,  100000  and  200000  cfs.   Inverse 
Gaussian  shape  hydrographs  with  ain  =  3,  6,  12  and  24  hours  were  routed 
through  a  rectangular  channel  (width  =  100  feet)  with  the  DWOPER 
model.   Transients  for  the  same  three  baseflows,  but  for  ain  "  6  and  12 
hours  were  routed  through  a  triangular  channel  with  side  slopes  of 
1:1.   The  bottom  slope  for  both  channels  was  1  foot  per  mile. 
Comparisons  between  the  numerical  solution  of  the  full 
Saint-Venant  equations  and  the  analytical  solution  of  equation  2.39  were 
made  at  100  and  400  miles  from  the  upstream  end  of  the  channels.  The 
values  of  c  and  D  in  equation  2.39  were  found  from  the  channel 
properties  at  each  baseflow  level.   Table  3.1  shows  the  combinations  of 
baseflow,  a*      and  channel  geometry  that  produce  the  results  for  small 
inflow  transients. 

It  should  be  emphasized  that  we  are  comparing  results  from  the 
numerical  solution  of  a  set  of  quasi-linear  partial  differential 
equations  with  results  from  the  analytical  solution  of  the  inverse 
Gaussian  pdf .   The  area  under  the  curve  defined  by  equation  2.39  will  be 
the  same  for  all  values  of  x.   The  area  under  a  hydrograph  is  the  volume 
of  water  in  the  transient.   Because  there  are  no  gains  or  losses  of 
water  in  the  prismatic  channels  studied,  the  volume  must  be  equal  for 
all  x.   This  is  guaranteed  by  equation  2.39. 
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Table  3.1   Organization  of  Figures  3.1  through  3.18, 
Comparing  the  Linearized  and  Full 
Saint-Venant  Equations 


Reference 

Flow 
(1000  cfs) 


10 


100 


200 


Rectangular  Channel 
ain  (hours) 

12 


24 


Figure  3.1  Figure  3.4  Figure  3.7  Figure  3.10 
Figure  3.2  Figure  3.5  Figure  3.8  Figure  3.11 
Figure  3.3   Figure  3.6   Figure  3.9   Figure  3.12 


10 

Reference   

Flow     100 
(1000  cfs)  


200 


Triangular  Channel 

°in  (hours) 
6      I      12 


Figure  3.13  Figure  3.16 
Figure  3.14  Figure  3.17 
Figure  3.15   Figure  3.18 
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Figures  3.1,  3.2  and  3.3  show  results  for  oin  =  3  hours  for 
baseflows  of  10000,  100000  and  200000  cfs,  respectively.   Outflows 
simulated  with  the  DWOPER  model  and  the  solution  of  the  inverse  Gaussian 
pdf  agree  closely  with  respect  to  the  timing  of  the  hydrographs  for  100 
and  400  miles  from  the  channel  inlet,  but  differ  in  the  magnitude  of  the 
flows.   In  these  cases  the  discrepancies  are  caused  by  problems  with  the 
numerical  solution  of  the  Saint-Venant  equations.   For  an  abrupt  wave, 
such  as  the  ain  =  3  hour  inflow  hydrograph,  the  distance  and  time  step 
sizes  for  the  numerical  solution  must  be  small.   For  the  results 
presented  in  Figures  3.1  through  3.3,  the  distance  step  was  2.5  miles 
and  the  time  step  was  0.25  hours.  With  so  many  distance  and  time  steps 
over  400  miles  and  100  hours,  numerical  roundoff  caused  the  DWOPER  model 
to  not  maintain  continuity.   Therefore  the  flows  predicted  by  the  DWOPER 
model  differ  from  those  found  with  equation  2.39.   However,  the  timing 
properties  of  the  very  abrupt  inflow  hydrograph  are  excellent. 

Slightly  less  severe  inflow  hydrographs,  a in  =  6  hours,  are 
depicted  in  Figures  3.4,  3.5  and  3.6.   Results  from  the  inverse  Gaussian 
pdf  and  the  DWOPER  model  compare  very  well,  especially  for  the  higher 
flows.   The  numerical  roundoff  problems  of  the  DWOPER  model  are  much 
less  evident  for  this  inflow  hydrograph. 

For  the  longer  duration  inflow,  with  oin  -  12  hours,  presented  in 

Figures  3.7,  3.8  and  3.9,  the  numerical  solution  of  the  full 

Saint-Venant  equations  is  matched  well  by  the  inverse  Gaussian  pdf.   For 

these  results,  as  well  as  those  shown  in  Figures  3.10,  3.11  and  3.12  for 

a        =24  hours,  the  analvtical  and  numerical  solutions  are  almost 
uin 

indistinguishable.  For  these  broader  hydrographs  the  distance  and  time 
step  sizes  can  be  increased  so  that  the  numerical  roundoff  is  no  longer 
a  problem. 
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To  demonstrate  the  ability  of  the  inverse  Gaussian  pdf  to  match 
the  DWOPER  model  for  a  general  prismatic  channel  shape,  defined  by 
equation  2.23,  results  for  ain  =  6  hours  in  a  triangular  channel  are 
presented  in  Figures  3.13,  3.14  and  3.15.   The  delay  in  the  arrival  of 
the  peak  flow  and  the  increased  dispersion  caused  by  the  triangular 
channel  shape  can  be  seen  by  comparing,  for  example,  Figure  3.4  with 
Figure  3.13.   Simply  by  changing  cQ  and  D  of  equation  2.39  to  account 
for  the  new  channel  shape,  the  inverse  Gaussian  pdf  remains  an  excellent 
approximation  to  the  DWOPER  model  when  the  inflow  transient  is  small. 

Figures  3.16,  3.17  and  3.18,  for  a ±n   =  12  hours  in  a  triangular 
channel,  can  be  compared  with  Figures  3.7  through  3.9.   Again  the  timing 
and  dispersive  properties  of  the  channel  are  identified  well  by  the 
inverse  Gaussian  pdf. 

Since  the  impulse  response  function  of  a  linear  system  is  a 
function  of  the  system  not  of  the  inputs,  cl.  and  D  have  the  same  values 
for  a  given  baseflow  for  all  rectangular  channel  results.   Similarly, 
for  the  triangular  channel,  cQ  and  D  are  functions  only  of  baseflow,  not 
°f  °in«   It  is  helpful  in  the  analysis  of,  for  example,  Figures  3.1, 
3.4,  3.7  and  3.10  (each  of  which  is  for  Q0  =  10000  cfs  in  a  rectangular 
channel),  to  realize  that  only  xQ  ,  as  defined  in  Figure  2.4,  is 
changing.   The  sets  of  results  for  the  same  baseflow  level  for  either 
the  rectangular  or  triangular  channel  shape  can  be  interpreted  as  an 
impulse  response  at  ever-increasing  distances  from  x  =  0,  the  location 
of  the  upstream  impulse. 

Figures  3.1  through  3.18  demonstrate  that,  for  a  wide  range  of 
flow  levels  and  hydrograph  durations  in  prismatic  channels  defined  by 
equation  2.23,  the  impulse  response  function  of  the  linearized 
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Saint-Venant  equations  (i.e.,  the  inverse  Gaussian  pdf)  produces  results 
almost  identical  to  those  from  the  DWOPER  model  for  small  inflow 
transients.   Since  results  from  the  inverse  Gaussian  pdf  are  found  by 
analytically  solving  equation  2.39,  the  CPU  time  required  is  small 
(i.e.,  less  than  a  second).   The  CPU  times  for  the  numerical  solution  of 
the  Saint-Venant  equations  on  a  PRIME  750  computer  are  presented  in 
Table  3.2.   The  DWOPER  model  used  from  1  to  3  orders  of  magnitude  more 
CPU  time  than  did  the  impulse  response  function  to  predict  the  flows 
shown  in  Figures  3.1  through  3.18.   In  the  next  section  the  impulse 
response  of  the  linear  state  space  model  is  compared  with  the  inverse 
Gaussian  pdf. 

Verification  of  the  Linear  State  Space  Model  for 
Small  Inflow  Transients 

Comparison  with  the  3-Parameter  Gamma  PDF 

The  three  parameters  of  the  linear  state  space  model  specified  by 
equations  2.91  through  2.98  are  a,  the  pure  delay;  n,  the  number  of 
reservoirs;  and  k,  the  time  delay  in  each  reservoir.   The  parameters  of 
the  3-parameter  gamma  pdf  (equation  2.69)  can  be  interpreted  identically 
to  the  parameters  of  the  state  space  model.   In  this  section,  results 
from  the  linear  state  space  model  and  the  3-parameter  gamma  pdf  are 
compared  when  the  same  values  of  a,  n  and  k  are  used  as  parameters  for 
both  the  state  space  model  and  the  pdf. 

Figure  3.19  depicts  the  strategy  used  to  verify  the  linear  state 
space  model.   The  parameters  c   and  D  of  the  inverse  Gaussian  pdf  can  be 
determined  for  a  specified  prismatic  channel  with  known  slope  and  shape 
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Table    3.2     CPU   Times   for   the   Numerical   Solution  of   the 

Saint-Venant   Equations   on  a  PRIME   750  Computer 
for   Small    Inflow  Transients 


10 


Reference      

Flow             100 
(1000   cfs)     


Rectangular  Channel 

ain   (hours) 

12 


24 


290.15  sec   156.65  sec    49.39  sec   20.90  sec 
215.23  sec   137.00  sec    40.50  sec   15.82  sec 


200  I   182.50  sec   127.08  sec    39.05  sec   15.75  sec 

■I 


10 

Reference   

Flow     100 
(1000  cfs)  


200 


Triangular  Channel 
oln  (hours) 

6      I      12 


144.56  sec  39.16  sec 
141.76  sec  30.74  sec 
128.87  sec    30.66  sec 
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Figure  3.19     Verification  Process  for  the 
Linear  State  Space  Model 


properties.   This  pdf  is  used  to  approximate  the  impulse  response  of  the 
Saint-Venant  equations  of  unsteady  flow.   Because  the  functional  form  of 
the  inverse  Gaussian  pdf  is  not  easily  expressed  in  state  space  form,  it 
is  approximated  by  a  3-parameter  gamma  pdf.   Parameters  of  the  two  pdfs 
are  related  by  matching  their  means  and  variances. 

The  degree  to  which  the  3-parameter  gamma  pdf  can  approximate  the 
inverse  Gaussian  pdf  has  been  shown  in  Figures  2.6  through  2.8.   Now  we 
demonstrate  that  the  linear  state  space  model  matches  the  gamma  pdf  when 
their  parameters  are  equal.   The  state  space  model  (like  any  numerical 
solution)  cannot  properly  simulate  given  an  impulse  input.  The  impulse 
must  be  averaged  over  a  computational  time  period,  At.   Although  At  can 
be  made  very  small,  that  increases  the  CPU  time  required  because  it 
involves  simulating  many  time  steps.   To  avoid  this  problem  while 
demonstrating  that  the  state  space  model  approximates  well  the 
3-parameter  gamma  pdf,  a  gamma  pdf  was  used  to  generate  a  smooth  inflow 
hydrograph.   Gamma  pdfs  can  be  added  if  the  parameter  k  is  equal  for  all 
the  distributions.   The  parameters  a  and  n  of  the  resultant  pdf  are 
found  by  summing  the  pure  delays  and  the  number  of  reservoirs, 
respectively,  of  the  gamma  pdfs  being  added.   The  parameters  of  the 
gamma  pdf  generating  the  inflow  are  defined  as  a± ,  r^  and  k-[ .   If  the 
parameters  of  the  linear  state  space  model  are  chosen  as  ara,  nm  and  km> 
the  expected  model  outflow  can  be  represented  as  a  gamma  pdf  with 

parameters  a0,  n0  and  k0  if 

a   =a.  +a  (3.3) 

o     l    m 

(3. A) 


n  =  n.  +  n 

o  i  ti 

k  =  k,  =  k 

o  i  n 


(3.5) 
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Equations  3.3  through  3.5  were  used  to  determine  parameters  for 
the  comparisons  shown  in  Figures  3.20  and  3.21.   The  same  a,  n  and  k 
values  used  in  Figures  2.6  and  2.7,  for  the  outflow  hydrograph,  are  used 
in  Figures  3.20  and  3.21,  respectively.  Values  for  the  parameters  a± 
and  ni  for  the  inflow  hydrograph  were  arbitrarily  chosen  as  0  and  2, 
respectively.   The  only  restrictions  on  ai  and  n-^  are 


and 


a4    <  a  (3-6) 

i  —  o 


n   <  n  (3.7) 

i    o 


The  close  agreement  between  the  linear  state  space  model  and  3-parameter 
gamma  pdf  shown  in  Figures  3.20  and  3.21  demonstrates  that  the  state 
space  model  is  an  excellent  approximation  of  the  3-parameter  gamma  pdf. 

Figures  2.6  and  2.7  show  how  well  the  inverse  Gaussian  and 
3-parameter  gamma  pdfs  match  when  their  parameters  are  related  with 
equations  2.73  through  2.75.   These  results,  and  the  excellent 
approximation  of  the  linear  state  space  model  for  the  3-parameter  gamma 
pdf  seen  in  Figures  3.20  and  3.21,  demonstrate  that  the  linear  state 
space  model  can  emulate  the  inverse  Gaussian  pdf. 

Comparison  with  the  DWOPER  Model 

The  next  step  in  Figure  3.19  is  to  compare  results  from  the  state 
space  and  DWOPER  models.   Small  (10  percent  of  baseflow)  transients  were 
routed  with  the  state  space  and  DWOPER  models. 

The  purpose  here  is  to  verify  that  outflow  hydrographs  from  the 
linear  state  space  model  compare  well  with  a  numerical  solution  of  the 
full  Saint-Venant  equations  when  the  inflow  transients  are  small.   To 
compare  these  two  models,  small  transients  were  input  to  the  DWOPER 
model.   The  mean  and  variance  were  computed  for  the  inflow  and  outflow 
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INFLOW  GAMMA  PDF 

DOTTED  LINE 

a  -  0,  n  -  2,  k  -  2.0 

OUTFLOW  GflMMfl  PDF 

DASHED  LINE  WITH  SYMBOLS 

a  -  1,  n  -  4,  k  -  2.0 

LINEAR  STATE  SPACE  MODEL 

SOLID  LINE 

a  "  I,   n  -  2,   k  -  2.0 


Figure  3.20   Comparison  of  the  Linear  State  Space 
Model  and  3-Parameter  Gamma  PDF 
for  a  =  1,  n  =  4  and  k  =  2.0 
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INFLOW  GAMMA  PDF 

DOTTED  LINE 

a  -  0,  n  -  2,  k   -  0.5 

OUTFLOW  GflMMfl  PDF 

DASHED  LINE  WITH  SYMBOLS 

a  -  8,   n  -  8,   k  -  0.5 

LINEAR  STATE  SPACE  MODEL 

SOLID  LINE 

a  -  8,   n  -  6,   k  -  0.5 


TIME 

Figure  3.21   Comparison  of  the  Linear  State  Space 
Model  and  3— Parameter  Gamma  PDF 
for  a  =  8,  n  =  8  and  k  =  0.5 
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hydrographs  of  the  DWOPER  model.   Parameters  a,  n  and  k  of  the 
3-parameter  gamma  pdf  (identical  to  the  parameters  of  the  linear  state 
space  model)  were  computed  from  these  two  moments  using  equations  2.78 
through  2.80.   These  three  parameters,  with  equations  2.91  through  2.98, 
define  the  linear  state  space  model.   Varying  the  value  of  p,  the  time 
weighting  factor,  from  0  to  1  had  no  impact  on  the  results  in  the  cases 
tested.   The  parameter  p  was  set  to  1  for  the  results  presented  to 
assure  stability  of  the  linear  equations  (Fread  1974,  pp.  7-10).   For 
purposes  of  this  comparison,  selected  channel  geometry  and  inflow 
combinations  from  the  previous  section  were  repeated.   Table  3.3  shows 
pairs  of  figures  with  identical  inflow  and  channel  properties. 

Now,  reimpose  the  integer  restriction  on  parameter  n.   In  the 
state  space  model,  n  represents  the  number  of  reservoirs  in  the  cascade 
and  must,  therefore,  be  a  positive  whole  number.   From  equations  2.70 
and  2.71,  the  mean  and  variance  of  the  3-parameter  gamma  pdf  (and  also 
of  the  impulse  response  function  of  the  state  space  model)  are 

U  =  a  +  n/k 
and 

a2  =  n/k2 
The  exact  values  of  a,  n  and  k  used  in  the  state  space  model  are 
presented  in  Table  3.4.   These  values  were  computed  from  the  moments  of 
the  inflow  and  outflow  hydrographs  of  the  DWOPER  model  using  equations 
2.78  through  2.80.   However,  n  must  be  an  integer.   The  values  of  n 
listed  in  Table  3.4  were  rounded  to  the  nearest  integer  for  use  in  the 
state  space  model. 

The  slight  effect  of  this  parameter  change  can  be  seen  by 
comparing  Figure  3.22  with  Figure  3.4.   The  ideal  value  of  n  is  3.702, 
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Table  3.3   Pairs  of  Figures  with  Identical  Input  Data 

Comparing  the  Linearized  Saint-Venant  Equations 
and  the  Linear  State  Space  Model  with  the  Full 
Saint-Venant  Equations 


Figure  Comparing  Figure  Comparing  Linear 

Linearized  and  Full  State  Space  Model  and  Full 

Saint-Venant  Equations  Saint-Venant  Equations 


3.4  3.22 

3.2  3.23 

3.8  3.24 

3.12  3.25 

3.14  3.26 

3.18  3.27 
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Table  3.4   Values  of  Parameters  a,  n  and  k 
for  the  State  Space  Model 


Rectangular  Channel,  L  =  100  miles 


10000  cfs 

a 

n 

k 

°^ 

5.692 

3.702 

0.3252 

°n 

= 

100000  cfs 

2.260 

0.943 

0.2087 

°o 

= 

200000  cfs 

1.711 

0.625 

0.1828 

Rectangular 

Channel, 

L  = 

400 

miles 

10000  cfs 

a 

n 

k 

On 

22.767 

14.810 

0.3253 

°n 

= 

100000  cfs 

9.039 

3.772 

0.2087 

% 

= 

200000  cfs 

6.844 

2.502 

0.1828 

Triangular 

Channel, 

L  = 

100 

miles 

10000  cfs 

a 

n 

k 

% 

6.830 

2.686 

0.1967 

% 

= 

100000  cfs 

3.841 

1.133 

0.1476 

% 

= 

200000  cfs 

3.230 

0.874 

0.1354 

Triangular 

Channel, 

L  = 

400 

miles 

10000  cfs 

a 

n 

k 

% 

27.318 

10.744 

0.1966 

% 

= 

100000  cfs 

15.363 

4.534 

0.1476 

% 

= 

200000  cfs 

12.919 

3.497 

0.1354 
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but  this  must  be  rounded  to  4.   Equations  2.70  and  2.71  indicate  that 
this  change  should  increase  the  lag  and  the  dispersion  introduced  by  the 
state  space  model.   This  is  exactly  what  occurs.  Whereas  in  Figure  3.4, 
the  timing  of  the  peak  for  the  impulse  response  was  identical  with  the 
DWOPER  model  and  the  peak  flow  slightly  high  at  100  miles,  the  peak  flow 
from  the  state  space  model  is  later  and  lower  in  Figure  3.22. 

The  rounding  of  parameter  n  can  sometimes  seem  beneficial,  as 
Figure  3.23  shows.   The  ideal  value  of  n  is  rounded  from  0.943  to  1. 
This  will  increase  the  lag  and  the  variance  of  the  state  space  model 
outflow  hydrograph.   By  slightly  delaying  the  rising  and  the  falling 
limbs  of  the  hydrograph,  the  state  space  model  flow  agrees  better 
visually  with  the  DWOPER  model  than  did  the  inverse  Gaussian  pdf  in 
Figure  3.2. 

A  more  expected  result  can  be  seen  by  comparing  Figures  3.24  and 
3.8.   In  this  case,  the  increased  variance  has  caused  the  peak  flow  to 
decrease  because  the  volume  of  water  in  the  hydrograph  is  fixed,  and  the 
falling  limb  of  the  hydrograph  for  the  state  space  model  has  been 

delayed. 

For  a  baseflow  of  200000  cfs  in  a  rectangular  channel,  n  would 
ideally  be  equal  to  0.625.   The  change  in  the  outflow  hydrograph  from 
the  impulse  response  function  of  Figure  3.12  to  the  state  space  model 
outflow  shown  in  Figure  3.25  is  caused  by  the  substantial  correction  of 
parameter  n  to  an  integer  value.   Since  the  maximum  time  shown  on  Figure 
3.25  is  larger  than  for  the  other  figures,  horizontal  differences 
between  curves  are  correspondingly  larger. 

A  case  of  n  being  rounded  down  to  the  next  lower  integer  can  be 
seen  in  Figure  3.26  for  a  triangular  channel.   In  this  case  the  lag  and 
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variance  of  the  state  space  model  decrease  because  n  is  less  than  the 
ideal.   Comparison  of  Figure  3.26  with  Figure  3.13  demonstrates  that  the 
state  space  model  outflow  is  distinctly  earlier  and  less  dispersed  than 
results  from  the  impulse  response  function.   This  causes  the  peak  flow 
for  the  state  space  model  in  Figure  3.26  to  be  higher  than  for  the 
DWOPER  model. 

State  space  model  results  for  a  baseflow  of  200000  cfs  in  a 
triangular  channel  are  presented  in  Figure  3.27.   The  effects  of 
increasing  n  from  0.874  to  1  can  clearly  be  seen  by  comparing  the 
outflow  hydrograph  with  Figure  3.18.   The  peak  flow  of  the  state  space 
model  in  Figure  3.27  is  lower  because  the  outflow  hydrograph  is  more 
dispersed.   The  increased  lag  can  be  seen  in  the  hydrograph's  falling 

limb. 

Results  from  the  linear  state  space  model  compare  very  well  with 

those  from  the  DWOPER  model.   Small  differences  in  results  can  be 
attributed  to  the  interger  restriction  on  parameter  n.   Differences  in 
CPU  times  required  to  simulate  with  the  two  models  are  large.   The  CPU 
times  on  a  PRIME  750  computer  for  the  six  figures  discussed  above  are 
presented  in  Table  3.5.   The  DWOPER  model  typically  took  from  2  to  3 
orders  of  magnitude  more  CPU  time  than  the  linear  state  space  model. 

The  Importance  of  Nonlinearities 

The  previous  results  show  that  the  linear  theory  developed  in 
Chapter  II  works  very  well  when  the  phenomena  being  modelled  are 
linear.   However,  in  many  (if  not  most)  of  the  surface  runoff  or  channel 
routing  problems  encountered,  the  input  fluctuations  are  large,  and 
nonlinear  behavior  is  observed.   In  this  section  the  linear  state  space 
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Table  3.5  CPU  Times  for  the  Linear  State  Space  and 
DWOPER  Models  on  a  PRIME  750  Computer  for 
Results  of  Figures  3.22  through  3.27 


Linear  State 


°in 

Length 

Channel 

% 

Space  Model 

DWOPER  Model 

Figure 
3.22 

(hours) 

6 

(miles) 
100 

Shape 
Rect. 

(cfs) 
10000 

CPU  Time 

CPU  Time 

0.21  sec 

161.57  sec 

3.23 

3 

100 

Rect. 

100000 

0.19  sec 

78.88  sec 

3.24 

12 

100 

Rect. 

100000 

0.09  sec 

45.63  sec 

3.25 

24 

100 

Rect. 

200000 

0.08  sec 

10.21  sec 

3.26 

6 

100 

Tri. 

100000 

0.13  sec 

47.32  sec 

3.27 

12 

100 

Tri. 

200000 

0.07  sec 

10.99  sec 
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model  is  compared  with  the  DWOPER  model  for  large  (2000  percent  of 
baseflow)  inflow  transients.   Since  the  linear  model  is  only  strictly 
applicable  for  a  small  flow  range  around  the  reference  discharge,  we 
expect  the  results  of  this  section  to  deterioriate  significantly  from 
those  presented  above.   If  this  were  not  the  case,  Chapter  IV  would  be 
titled  "Conclusions"  instead  of  "A  Nonlinear  State  Space  Model." 

The  cases  presented  are  for  inflow  hydrographs  of  200000  cfs  on  a 
baseflow  of  10000  cfs.   Results  from  the  linear  state  space  and  DWOPER 
models  are  compared  at  100  and  400  miles  from  the  channel  inlet.   The  16 
cases  presented  are  summarized  in  Table  3.6. 

As  shown  in  Figures  3.28  and  3.29,  a  linear  model  is  not  adequate 
to  simulate  the  physical  processes  that  occur  when  a  large,  very  abrupt 
hydrograph  (oin  =  3  hours)  is  input.   Since  results  from  the  linear 
state  space  model  are  sensitive  to  the  reference  flow  level  at  which  the 
parameters  a,  n  and  k  are  selected,  Q0  was  chosen  as  100000  cfs  for  all 
large  transient  figures  unless  otherwise  stated.   This  value  was  chosen 
at  slightly  less  than  the  average  of  the  base  and  peak  flow  levels,  to 
assure  that  values  for  the  parameters  of  the  state  space  model  are  as 
often  too  high  as  too  low.   The  results  of  Figure  3.28  indicate  that, 
even  at  100  miles  from  the  channel  inlet,  the  state  space  model  does  not 
disperse  the  inflow  transient  as  much  as  the  DWOPER  model.   At  400 
miles,  shown  in  Figure  3.29,  outflow  from  the  state  space  model  is 
clearly  late  and  not  as  dispersed  as  from  the  DWOPER  model. 

Figure  3.30  shows  that  the  timing  and  dispersion  properties  of  the 
state  space  model  for  oln  =  6  hours  at  L  =  100  miles  are  very  close  to 
those  of  the  DWOPER  model;  however,  the  rising  and  falling  limbs  of  the 
state  space  model  outflow  hydrograph  have  the  wrong  slopes.   This  is  a 
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Table  3.6   Organization  of  Figures  3.28  through  3.43, 
Comparing  the  Linear  State  Space  Model 
with  the  Full  Saint-Venant  Equations 
for  Large  Inflow  Transients 


Rectangular  Channel 


'in 
'in 


=  3  hours 


'in 


'in 


12  hours 


L  =  100  miles 

Figure  3.28 

Figure  3.30 

Figure  3.32 

Figure  3.34 


L  =  400  miles 

Figure  3.29 

Figure  3.31 

Figure  3.33 

Figure  3.35 


Triangular  Channel 


a,      -      6  hours 
uin 

°"in  =  12  hours 


L  =  100  miles 

Figure  3.36 
Figure  3.38 


L  =  400  miles 

Figure  3.37 
Figure  3.39 


Double-Peaked  Hydrograph 


00  =  100000  cfs 
Qq  =  10000  cfs 
0Q  =  200000  cfs 


L  =  100  miles 
Figure  3.40 


>   =  400  miles 

Figure  3.41 
Figure  3.42 
Figure  3.43 
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typical  problem  with  linear  models  trying  to  simulate  nonlinear 
behavior.   The  linear  results  are  closer  to  a  Gaussian  pdf  than  to  the 
truely  nonlinear  shape  of  the  DWOPER  model  outflow  hydrograph.   This 
will  be  a  recurring  problem  for  the  linear  state  space  model  results 
presented  in  the  remainder  of  this  chapter. 

Results  for  a ±n   =  6  hours  at  L  -  400  miles,  shown  in  Figure  3.31, 
demonstrate  that  it  was  merely  good  fortune,  not  a  properly  formulated 
model,  that  let  the  linear  state  space  model  match  the  DWOPER  model 
outflow  in  Figure  3.30.   At  400  miles,  the  timing  and  dispersive 
properties  of  the  state  space  model  are  clearly  different  than  the 
nonlinear  results  from  the  DWOPER  model. 

With  the  broader  inflow  hydrograph  (oin  =  12  hours)  shown  in 
Figure  3.32,  the  behavior  of  the  state  space  model  improves  greatly. 
Since  the  flow  varies  more  gradually  over  time,  the  nonlinear  terms  of 
the  Saint-Venant  equations  are  relatively  less  important  than  for  the 
more  rapidly  varying  inflows  shown  earlier.   The  major  discrepancies 
between  the  state  space  and  DWOPER  models  are  the  slight  increase  in  lag 
and  dispersion  of  the  inflow  hydrograph  shown  by  the  state  space  model 
which  is  caused  by  rounding  parameter  n  up  to  1,  and  the  symmetrical 
shape  of  the  state  space  model's  outflow  hydrograph,  which  is  typical  of 

linear  models. 

The  Gaussian  shape  of  the  state  space  model  outflow  is  even  more 
evident  in  Figure  3.33,  for  a±n   =  12  hours  and  L  =  400  miles.   For  this 
slowlv  varying  inflow,  the  delay  and  dispersion  introduced  by  the  state 
space  model  are  approximately  equal  to  those  seen  in  the  DWOPER  model 
results,  but  the  linear  model  is  not  able  to  reproduce  the  slopes  of  the 
steeply  rising  and  slowly  falling  outflow  hydrograph  limbs  which  are 
characteristic  of  nonlinear  routing  models. 
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As  the  dispersion  of  the  inflow  hydrograph  increases,  these 
effects  become  even  more  pronounced.   In  Figures  3.34  and  3.35,  for 
L  =  100  and  400  miles,  respectively,  the  basic  timing  and  attenuation 
properties  of  the  state  space  and  DWOPER  models  are  in  close  agreement, 
but  the  symmetrical  shape,  typcial  of  the  linear  model,  limits  its 
ability  to  match  the  nonlinear  outflow  hydrograph. 

For  a  triangular  channel,  as  shown  in  Figure  3.36,  the  major  cause 
of  differences  between  the  state  space  and  DWOPER  models  is  the  rounding 
of  parameter  n  down  to  1.   As  we  will  see  in  Chapter  IV,  the  correct 
value  for  m,  the  kinematic  wave  parameter,  is  4/3  for  a  triangular 
channel.   This  is  close  to  m  =  1,  the  value  implied  by  a  linear  model. 
Since  m  should  be  equal  to  5/3  for  a  rectangular  channel,  we  would 
expect  a  linear  model  to  be  better  for  a  triangular  channel  than  for  a 
rectangular  one.   This  is  indeed  the  case,  as  can  be  seen  by  comparing 
Figure  3.37,  for  a^n  =   6  hours  and  L  =  400  miles  in  a  triangular 
channel,  with  the  corresponding  Figure  3.31  in  a  rectangular  channel. 

The  shape  of  the  state  space  model  outflow  in  Figure  3.37  is 
closer  to  the  DWOPER  model  than  in  Figure  3.31.   This  is  not  because  the 
state  space  model  is  any  less  linear  in  Figure  3.37,  but  because  flow  in 
a  triangular  channel,  with  m  =  4/3,  is  more  linear  than  in  a  rectangular 
channel,  where  m  =  5/3.   This  effect  is  even  more  evident  in  Figures 
3.38  and  3.39  for  an  inflow  hydrograph  with  a ln  =  12  hours  at  L  =  100 
and  400  miles,  respectively.   The  state  space  model  outflow  hydrograph 
shapes,  especially  the  falling  limbs,  are  much  closer  to  the  nonlinear 
behavior  simulated  with  the  DWOPER  model  than  in  the  corresponding 
Figures  3.32  and  3.33,  for  a  rectangular  channel. 
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The  next  test  comparing  the  linear  state  space  and  DWOPER  models 
was  for  a  double-peaked  inflow  hydrograph.   The  inflow  hydrograph 
consisted  of  a  100000  cfs  transient  followed  by  a  200000  cfs  transient 
on  a  baseflow  of  10000  cfs.  At  100  miles  from  the  inlet  of  a 
rectangular  channel,  the  linear  state  space  model  still  approximates  the 
DWOPER  model  results  well.   However,  as  shown  in  Figure  3.40,  the  rising 
and  falling  linbs  of  the  linear  model's  hydrograph  are  already  different 
in  slope  than  the  DWOPER  model.   In  Figure  3.41  we  can  see  that,  by  400 
miles  from  the  channel  inlet,  the  state  space  model  outflow  has  departed 
significantly  from  the  DWOPER  model.   The  outflow  hydrograph  from  the 
linear  model  is  becoming  more  symmetrical,  while  the  rising  limb  of  the 
DWOPER  model  outflow  is  steepening  and  the  falling  limb  is  becoming 
flatter. 

The  impact  of  the  reference  flow  level  on  the  linear  state  space 
model  results  is  shown  in  Figures  3.42  and  3.43.   The  parameters  a,  n 
and  k  were  selected  at  a  reference  flow  of  10000  cfs  in  Figure  3.42. 
The  resultant  outflow  hydrograph  is  very  late  and  much  more  dispersed 
than  the  flows  shown  in  Figure  3.41.   When  parameters  are  chosen  at 
q  =  200000  cfs,  as  in  Figure  3.43,  the  state  space  model  lags  and 
disperses  the  inflow  hydrograph  significantly  less  than  when  0Q  =  100000 
cfs.   The  reference  flow  level  has  a  striking  effect  on  the  linear  state 
space  model  behavior. 

Summary 

Results  presented  in  this  chapter  have  demonstrated 
1.   the  behavior  of  the  impulse  response  function  for  the 

linearized  Saint-Venant  equations  (equation  2.39)  is  almost 
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identical  with  a  numerical  solution  of  the  full  Saint-Venant 
equations  for  small  (10  percent  of  baseflow)  inflow 
transients , 

2.  the  linear  state  space  model  developed  in  Chapter  II  is  an 
excellent  surrogate  for  a  model  based  on  the  full 
Saint-Venant  equations  when  small  inflow  transients  are 
routed,  and 

3.  important  nonlinearities  that  are  not  adequately  modelled  by 
the  linear  theory  developed  in  Chapter  II  when  large  (2000 
percent  of  baseflow)  inflow  hydrographs  are  routed  with  the 
full  Saint-Venant  equations. 

In  Chapter  IV  nonlinearities  important  when  routing  large  inflow 
hydrographs  are  analyzed.   A  nonlinear  state  space  model  emulate  the 
behavior  of  the  full  Saint-Venant  equations  over  large  flow  ranges  is 
developed.   The  ability  of  this  nonlinear  state  space  model  to  act  as  a 
surrogate  for  the  DWOPER  model  is  the  topic  of  Chapter  V. 


CHAPTER  IV 
A  NONLINEAR  STATE  SPACE  ROUTING  MODEL 


Analysis  of  the  Important  Nonlinearltles 

Chapter  III  presented  limitations  of  the  linear  theory.   The 
inability  of  the  linear  theory  (specifically  the  impulse  response 
function  of  the  linearized  Saint-Venant  equations)  to  emulate  the 
behavior  of  the  full  nonlinear  Saint-Venant  equations  leads  us  to  the 
current  chapter.   In  this  chapter  we  will 

1.  analyze  the  nonlinear  characteristics  of  the  primary  model, 

2.  explore  the  properties  of  a  nonlinear  storage  reservoir,  and 

3.  develop  a  fully  nonlinear  state  space  routing  model  which 
emulates  the  mathematical  behavior  of  the  Saint-Venant 
equations  for  any  inflow  transient. 

The  approach  to  the  analysis  of  the  nonlinearities  inherent  in  the 
Saint-Venant  equations  is  outlined  in  Figure  4.1.   The  following  section 
begins  with  the  full  Saint-Venant  equations  and,  in  the  spirit  of  the 
new  perspective  of  routing  presented  in  Figure  1.1,  predicts  the 
theoretical  variation  of  the  impulse  response  properties  with  a  linear 
analysis.   The  expectations  of  linear  theory  are  then  confirmed  by 
numerically  solving  the  fully  nonlinear  Saint-Venant  equations  and 
observing  the  changes  in  the  impulse  response  function.   Results  from 
this  analysis  lend  insight  into  a  state  space  model  structure  which  will 
emulate  the  nonlinear  behavior  of  the  primary  model. 
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Figure  4.1    Approach  to  the  Nonlinear  Analysis 


125 


Nonlinear  Characteristics  of  the  Primary  Model 


Introduction 


This  section  analyzes  the  nonlinear  characteristics  of  the 
Saint-Venant  equations,  the  exemplar  primary  model.   The  approach  for 
that  analysis  is  shown  in  Figure  4.2. 

The  variations  with  flow  of  the  linearized  Saint-Venant  equations' 
impulse  response  function  are  derived  and  the  mean  and  variance 
properties  of  the  impulse  response  function  are  analyzed.   These  moments 
will  be  needed  to  specify  the  parameters  of  the  nonlinear  state  space 
model  developed  later. 

The  variations  with  flow  of  the  impulse  response  for  the  full 
Saint-Venant  equations  are  determined  by  simulating  with  incremental 
transients  as  inflow,  and  observing  the  outflow.   The  results  predicted 
by  the  linear  analysis  of  the  Saint-Venant  equations  agree  with  those 
obtained  by  simulating  with  the  DWOPER  model. 

Theoretical  Variations  with  Flow  Level  of  the  Impulse  Response 
Function  for  the  Linearized  Saint-Venant  Equations  for  Incremental 
Inflow  Transients 

Beginning  with  the  Saint-Venant  equations  of  unsteady  flow 
(equations  2.20  and  2.21),  the  impulse  response  function  for  the 
linearized  Saint-Venant  equations  was  obtained  in  Chapter  II  as 
equation  2.39 

(x-c  t)2 

4'TT«D't3 
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This  equation  has  the  form  of  an  inverse  Gaussian  pdf  (Johnson  and  Kotz 
1970,  p.  137).  The  mean,  t^,  and  variance,  a2,  of  the  inverse  Gaussian 
pdf  were  presented  in  equations  2.41  and  2.42  as 


and 


'i-f 


a2  =  2«  —  fi-V 


S  -c3 

0 

The  derivation  of  the  changes  with  flow  level,  0,  in  the  mean  and 
variance  of  the  linearized  Saint-Venant  equations'  impulse  response 
function  begins  by  substituting  for  c  in  equation  2.41  with  equation 
2.38  to  give 


t.   --*-  (4.1) 

x.    ra*  v 
From  equation  2.22 


so 


t.  =-•-  (4-2) 

I        m  0 


Equation  2.30  expresses  the  friction  slope  as 


2 

S  =_^ ,S 

f    2  ,  2m  n 
a  'A    u 


Low 


Since  this  analysis  is  for  small  transients  about  a  reference  flc 
level,  we  make  the  standard  kinematic  assumption  of  steady  uniform  flow 

Sc  =  S  (4.3) 

f     0 

(Henderson  1966,  pp.  287,367).   Substituting  equation  4.3  into  equation 
2.30  and  rearranging  gives 


(i)1/m-Q1/m 


Replacing  A  in  equation  4.2  with  equation  4.4  produces 
1/m 


sx,fM/B,Q 
A        m  ^^  0 
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(4.4) 

(4.5) 


1-m 


h  =  vQ 


where 


x  f  U  1/m 


Recognizing  that  8.  is  constant  for  a  reference  flow 
1-m 

tt~   0 


(4.6) 


(4.7) 


where  ~  symbolizes  "is  proportional  to." 

Equation  4.7  states  that  the  first  moment  of  the  impulse  response 
function  varies  as  a  simple  algebraic  power  of  the  flow  level.   The 
value  of  the  exponent  depends  upon  the  kinematic  wave  parameter,  m. 

The  variation  of  the  second  moment  can  be  found  by  again 
substituting  for  c  with  equation  2.38.   This  time  substitution  into 
equation  2.42  gives 


,2  „ 


Q-x 


B«  S  •  m3  •  v3 

0 


1H>' 


(4.8) 


where   0  f    bq  and  cQ  are  now  replaced  by  Q,  B  and  c  because  we  are 

solving  for  the  variation  of  a2  with  flow. 
Substituting  for  v  with  equation  2.22  produces 


B'S  «m3   Q3 


(4.9) 


129 


By  replacing  B  with  equation  2.23  we  get 

a2  =  9^. ,(*!).(  i-*2)  (4.10) 

B  «Ar»S  •m3   03 
0     0 

And,  finally,  substituting  for  A  with  equation  4.4  and  rearranging 
yields 

3-r    £2- 


t2  „ 


I   •  S  •  m3 
0   0 


-(iV)-Ci-)  m  '^—  (4.11) 


By  simplifying  the  expression  for  Q  we  obtain 

3-r- 2m 

a2  =  0  -Q   m  (4.12) 

2 

where 

m 

e    = l_-(i-*2)-ar-3 

2   B  •  S  •  ra3 

0   0 

The  standard  deviation  of  the  impulse  response  is 
3-r- 2m 

a   =/-r  »  /b--q  2m  (4.13) 


3-r-2m 


a  ~  Q 


2m  (4.14) 


Again,  the  moment  of  the  impulse  response  varies  as  a  simple  algebraic 
power  of  0  and  the  exponent  depends  on  m.   For  the  variance,  however, 
the  power  is  also  a  function  of  the  channel  shape  parameter,  r.   For  a 
given  friction  resistance  formula  (i.e.,  Manning  or  Chezy) ,  m  and  r  are 
related  by  equation  2.28. 

The  coefficient  of  variation  of  the  linearized  Saint-Venant 
equations'  impulse  response  can  be  found  from  equation  2.43 


130 


a 

Substituting  equations  4.6  and  4.13  into  equation  2.43  gives 

3-r-2m 
/3     2m 

%=^Q-T^  (4'15) 

1   Q~ 


or  simplifying 


1-r 

c  =  6  -Q  2m  (4.16) 

v    3 

where 

/3~ 

Q   _    2 

33    X 
1 

or  finally 

1-r 

c  ~  0  2m  (4.17) 

v 

Equations  4.7,  4.14  and  4.17  express  the  theoretical  variations  with 
flow  level  of  the  mean,  standard  deviation  and  coefficient  of  variation 
of  the  impulse  response  function  of  the  Saint-Venant  equations.   In  the 
following  section  some  numerical  results  are  presented  for  comparison 
with  the  variations  predicted  by  linear  theory. 

Simulated  Variations  with  Flow  Level  of  the  Impulse  Response  Function 
for  the  Full  Saint-Venant  Equations  for  Incremental  Inflow  Transients 

In  general,  the  operational  procedure  for  determining  the  moments 
of  the  channel  response  will  be  to  simulate  small  transients  with  the 
primary  model  at  various  baseflow  levels  and  observe  the  results.   In 
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the  current  application,  however,  with  the  Saint-Venant  equations  as  the 
primary  model  for  a  prismatic  channel,  we  have  the  opportunity  to 
compare  the  expected  theoretical  results  obtained  from  equations  4.7, 
4.14  and  4.17,  with  numerical  computation  of  the  moments  from  the  fully 
nonlinear  Saint-Venant  equations  model. 

The  NWS  DWOPER  model  (Fread  1978)  was  used  to  numerically  solve 
the  full  equations  of  unsteady  flow.   The  moments  of  the  impulse 
response  function  of  the  Saint-Venant  equations  were  found  by  driving 
the  DWOPER  model  with  small  (10  percent  of  baseflow)  inflow  transients 
for  several  baseflow  levels,  and  observing  the  outflows.   The  means  and 
variances  of  the  inflow,  x,  and  outflow,  y,  hydrographs  were  then 
computed.   The  moments  of  the  impulse  response  function,  h,  were  found 
with  equations  2.53  and  2.54 


U'(h)  =  u'(y)  -  U'(x) 


and 


U  (h)  =  U  (y)  -  U  (x) 

2       2       2 


The  coefficient  of  variation  of  the  impulse  response  is 


/U  (h) 
c   =   2,  (4.18) 

Cv    U'(h) 

1 

From  equations  4.7,  4.14  and  4.17,  we  see  that  the  theoretical 
variations  with  flow  of  the  impulse  response  moments  depend  on  m  and 
r.   The  theoretical  values  of  m  can  be  derived  for  any  prismatic  channel 
using  the  Manning  friction  formula.   Once  m  is  known,  the  value  of  r  can 
be  found  from  equation  2.28.   For  the  two  channel  shapes  studied  in  this 
work,  the  values  of  m  and  r  are: 
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shape  m  r 

rectangular  5/3  0 

triangular  4/3  1/2 

Figures  4.3  and  4.4  depict  the  lag  time  introduced  by  the 
rectangular  and  triangular  channels,  respectively.   The  geometry  of  the 
channels  is  described  in  Chapter  III.   The  lags  presented  are  for 
distances  of  25,  50,  100,  200  and  400  miles  from  the  channel  inlet  (Xq 
in  the  notation  of  Chapter  III).   Incremental  transients  on  baseflows  of 
10000,  100000  and  200000  cfs  were  routed  with  the  DWOPER  model.   The 
results  presented  are  for  inflow  transients  with  o"in  =  12  hours. 
Results  for  a..   ■  3,  6  and  24  hours  look  identical  to  those  presented. 
This  must  be,  because  the  impulse  response  for  a  linear  system  is  not 
dependent  on  the  properties  of  the  input,  but  is  strictly  a  function  of 
the  system.   The  theoretically  predicted  and  numerically  simulated 
values  match  almost  exactly.   The  slope  of  the  lines  on  Figures  4.3  and 
4.4  can  be  found  by  rearranging  equation  4.7  as 

slope,    -  ™  (4.19) 

lag   1-m 

For  the  rectangular  channel  presented  in  Figure  4.3,  with  m  =  5/3,  the 

slopei    =  -2.5.   On  Figure  4.4,  the  triangular  channel  with  m  =  4/3, 

the  slope^   =  -4.   Observed  lag  time  versus  flow  data  which  corroborate 

the  theoretical  expectations  of  equation  4.19  are  presented  in  Chapter 

VII. 

Figures  4.5  and  4.6  present  the  theoretical  and  simulated  standard 

deviations  of  the  impulse  response  function  versus  flow  level.   Again, 

the  results  compare  very  well.   The  slopes  of  the  lines  on  Figures  4.5 

and  4.6  can  be  found  by  substituting  (m  =  5/3,  r  =  0)  and  (m  =  4/3, 

r  =  1/2),  respectively,  into 
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(4.20) 


Similarly,  Figures  4.7  and  4.8  depict  the  theoretical  and  simulated 
coefficients  of  variation  of  the  impulse  response  function.   The  slopes 
of  these  lines  can  be  found  as 


slope 


2m 
1-r 


(4.21) 


Before  leaving  the  topic  of  the  simulated  variance  with  flow  level  of 
the  impulse  response  function,  a  valuable  relationship  can  be  found  by 
solving  equation  4.19  for  m  as 


slope 


lag 


1  +  slope.. 

*   lag 


(4.22) 


This  expression,  defining  the  kinematic  wave  parameter,  ra,  in  terms  of 
the  channel  lag  versus  flow  properties,  will  be  useful  for  determining 
one  of  the  parameters  of  the  nonlinear  state  space  model  which  is 
developed  in  a  subsequent  section  of  this  chapter. 


Properties  of  the  Nonlinear  Storage  Reservoir 
as  a  Building  Block  for  the  State  Space  Model 

Just  as  the  linear  reservoir  relationship,  equation  2.82,  is  the 
cornerstone  of  the  linear  state  space  model,  the  nonlinear  reservoir 
equation  is  the  foundation  in  the  development  of  the  nonlinear  state 
space  routing  model.   The  relationship  between  storage  and  outflow  for  a 
nonlinear  reservoir  is 


S   =  <•  0 


where   S  =  the  reservoir  storage, 


(4.23) 
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0  =  the  reservoir  outflow,  and 

m  and  <  are  constants. 
If  m  in  equation  4.23  is  taken  to  be  the  kinematic  wave  parameter, 
the  nonlinear  storage  reservoir  has  the  same  lag  versus  flow  properties 
as  the  linearized  Saint-Venant  equations.   Equation  4.23  can  be 
rewritten  as 

s  =  <1/m.o1/m  (4.24) 

For  a  reservoir  of  length  L  and  cross  sectional  area  A 

S  =  L-A  (4.25) 

Combining  equations  4.24  and  4.25  gives 

,    1   1/m   1/m  ,.  ... 

A  =  j-k        «0  (4.26) 

The  time  lag  through  a  reservoir  is 

t.  =1  (4.27) 

Jc   v 

where   v  is  the  average  velocity  through  the  reservoir. 
But  from  equation  2.22 

0 

V  "  A 

Substituting  equation  2.22  into  equation  4.27  produces 


h    =  L4  (4'28) 

And  substituting  equation  4.26  into  equation  4.28  yields 
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1-m 
t  =  <1/m.0  m  (4.30) 

which  is  identical  in  form  to  equation  4.6  for  the  time  lag  of  the 
linearized  Saint-Venant  equations. 

The  thesis  presented  in  this  work  is  that  a  state  space  model  in 
the  form  of  a  cascade  of  nonlinear  storage  reservoirs  can  emulate  the 
nonlinear  behavior  of  the  Saint-Venant  equations.  A  fully  nonlinear 
state  space  model  that  maintains  the  lag  versus  flow  relationship  of 
equations  4.6  and  4.30  is  developed  in  the  next  section. 


Development  of  the  Nonlinear  Model 


Introduction 


A  fully  nonlinear  state  space  routing  model  which  emulates  the 
Saint-Venant  equations  for  downstream  waves  in  prismatic  channels  is 
developed  in  this  section.   This  nonlinear  state  space  model  is  based 
upon 

1.  the  linear  state  space  model  developed  in  Chapter  II,  which 
has  the  same  Impulse  response  function  as  the  linearized 
Saint-Venant  equations  for  a  prismatic  channel  (i.e.,  the 
Inverse  Gaussian  pdf);  and 

2.  the  functional  relationships  for  the  variations  with  flow  of 
the  mean  and  variance  of  the  inverse  Gaussian  pdf  which  were 
derived  earlier  in  this  chapter. 
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Figure  4.9  shows  the  approach  taken  in  developing  the  nonlinear 
state  space  model.   The  general  structure  of  the  nonlinear  state  space 
model  is  a  cascade  of  nonlinear  reservoirs  as  shown  in  Figure  4.10.   The 
nonlinear  properties  are  implemented  by  varying,  with  flow  level,  the 
parameters  (a,  n  and  k)  of  the  linear  state  space  model  developed  in 
Chapter  II.   Instead  of  a  model  specified  by  three  parameters,  as  In 
Chapter  II,  the  nonlinear  state  space  model  is  specified  by  three 
functions.   The  values  of  a,  n  and  k  vary  with  the  level  of  flow. 
Although  each  of  the  n  reservoirs  in  the  cascade  is  governed  by 
identical  functions  of  a  versus  flow  and  k  versus  flow,  in  general  the 
values  of  a-i  and  k.<  will  be  different  for  each  ith  reservoir  because, 
except  under  steady  flow  conditions,  the  flow  through  each  reservoir 
will  be  different.   From  three  functional  relationships,  2n+l  parameters 
are  determined  which  specify  the  F^   t  and  _Gj(  t  matrices  in  equation  1.1 
at  each  time  step.   The  algorithm  to  perform  this  task  efficiently  is 
presented  later  in  this  chapter,  but  first  the  functions  of  a,  n  and  k 
versus  flow  are  derived. 

Equations  4.7  and  4.14  express  the  variations  with  flow  of  the  lag 
and  standard  deviation  of  the  impulse  response  function  for  the 
linearized  Saint-Venant  equations.   We  want  the  variations  with  flow  of 
the  lag  and  variance  properties  of  the  state  space  model  to  match  the 
variations  of  the  fully  nonlinear  Saint-Venant  equations. 

The  parameters  of  the  state  space  model  are  related  to  the  mean 
and  variance  of  the  impulse  response  function  by  equations  2.78  through 
2.80 
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ASSUME  STRUCTURE  OF  STRTE  SPRCE  MODEL 

EQUIVALENT  TO  fl  CASCADE  OF  NONLINEAR 

RESERVOIRS  PLUS  A  PURE  DELAY  THAT 

VARIES  WITH  FLOW  LEVEL 


DERIVE  ANALYTICAL  EXPRESSIONS  FOR 

MOMENTS  OF  THE  IMPULSE  RESPONSE 

FUNCTION  OF  THE  STRTE  SPACE  MODEL 

AS  A  FUNCTION  OF  FLOW  LEVEL 


i 


EXPRESS  PARAMETERS  OF  3-PARRMETER 

GAMMA  PDF  AS  A  FUNCTION  OF  MOMENTS 

AND  FLOW  LEVEL 
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EXPRESS  STATE  TRANSITION  AND  INPUT 

COEFFICIENT  MATRICES  AS  FUNCTIONS 

OF  FLOW  LEVEL 


Figure  4.9    Approach  to  the  Development  of  the 
Nonlinear  State  Space  Model 
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Figure  4.10  A  Cascade  of  Nonlinear  Reservoirs 
Governed  by  S111  =  k*0 
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1 
a  =T* 

PIG 

4 
n  =  9   ' 

°IG 

-!• 

a2 
IG 

Substituting  equation  4.7  into  equation  2.78  with  t   =  u    gives  the 

Jo     IG 

variation  of  the  state  space  parameter  a  with  flow  as 

1-m 

a~  Q  m  (4.31) 

The  variation  of  n  with  flow  is  found  by  substituting  equations  4.7  and 

4.14  into  equation  2.79  with  t.  =  u   and  a  =  a    to  give 

Ji     LG  1G 

2-2m 


-2 «-  (4.32) 

3-r-2m 


or 

r-1 
n  ~  0  m  (4.33) 

Similarly,  the  variation  of  k.  with  flow  can  be  obtained  from  equations 

4.7,  4.14  and  2.80  as 

1-m 

m 

k~   ,Q   ,,  (4.34) 

j-r-2m 

m 
0 


or 


m-2+r 
k~  0  m  (4.35) 
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Equations  4.31,  4.33  and  4.35  describe  the  variations  with  flow  of 
parameters  a,  n  and  k,  respectively,  of  the  state  space  model. 

Operationally,  the  value  of  m  for  the  channel  can  be  found  by 
simulation.   Using  the  primary  model,  the  lags  for  several  baseflows 
with  small  input  transients  can  be  found.   The  slope  of  the  line 
obtained  when  lag  is  plotted  as  a  function  of  baseflow  on  logarithmic 
axes  determines  m  with  equation  4.22. 

The  lag  and  variance  produced  by  the  channel  must  be  found  for  any 
one  flow  level  by  simulating  with  the  primary  model.  With  this  as  a 
reference  point  in  lag-Q  space,  ra  computed  as  above,  and  the  functions 
4.31,  4.33  and  4.35  for  the  model  parameters,  all  the  information  needed 
to  construct  the  nonlinear  state  space  model  is  at  hand.   The  steps 
required  to  develop  the  nonlinear  state  space  model  are  summarized  in 
Figure  4.11. 

Two  points  of  a  general  nature  should  be  stressed  at  this  time 

1.  The  state  space  model  developed  here  is  fully  nonlinear. 
Its  theoretical  basis  draws  heavily  from  linear  system 
analysis,  but  only  to  explain  better  the  nature  of  the 
important  nonlinearities;  and 

2.  There  is  not  a  reference  discharge  for  this  model.   The 
functions  described  by  equations  4.31,  4.33  and  4.35  are 
valid  for  all  flows  where  m  is  constant.   A  point  in  lag-Q 
space  to  fix  these  functions  could  be  found  at  any  flow 
level. 

The  following  sections  describe  the  final  sucessful  nonlinear 
state  space  model  and  the  path  taken  to  find  the  final  model.   The  first 
is  needed  to  confirm  the  perspective  of  this  work.   The  second  is 


147 


RESERVOIR 

CONTINUITY 

EQUATION 

AS/Ai  -  I  -  0 


NONLINEAR 

RESERVOIR 

EQUATION 

S-  -  c-0 


GENERAL  STATE 
SPACE  EQUATION 


CHANNEL  IMPULSE 

RESPONSE  FUNCTION 

AT  SEVERAL 

REFERENCE  FLOW 

LEVELS 


DETERMINE  FLOW  RANGES  IN 

WHICH  m  IS  CONSTANT  AND 

COMPUTE  a0,  n0  AND  k0  AT  A 

REFERENCE  FLOW  IN  EACH  RANGE 


EQUATIONS  FOR 

NONLINEAR 

RESERVOIRS 

IN  A  CASCADE 


i 


NONLINEAR  STATE  SPACE 
MODEL  WITH  FX1  AND  Gxt 
AS  FUNCTIONS  <fr  CHANNEL 
IMPULSE  RESPONSE 


Figure  4.11    Development  of  the  Nonlinear 
State  Space  Model 
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included  because  it  may  lend  insight  into  the  alternatives  available 
even  after  the  theory  is  developed.   Obviously,  both  the  theoretician 
and  the  model  builder  are  vital  if  a  useful  product  is  to  be  developed. 

The  Final  Configuration  of  the  Nonlinear  Model 

The  final  form  of  the  nonlinear  state  space  routing  model  is  based 
conceptually  on  a  cascade  of  nonlinear  reservoirs.   As  with  the  linear 
state  space  model  developed  in  Chapter  II,  equation  2.81  describes  the 
relationships  among  inflow,  outflow  and  storage  for  a  reservoir 

dT=  I"° 

where 

S  =  the  reservoir  storage, 

I  =  the  reservoir  inflow,  and 

0  =  the  reservoir  outflow. 
Now,  however,  the  relationship  between  storage  and  outflow  is  nonlinear, 
as  expressed  in  equation  4.23.   The  nonlinearites  implicit  in  equation 
4.23  are  modelled  by  letting  k  vary  with  flow  for  the  ith  reservoir  as 

S    =  K   •()  ,  (4.36) 

t,i    t,i   t,i 

where   the  subscript  t,i  indicates  the  value  of  <    for  the  flow  level 
through  the  ith  reservoir  during  the  tth  time  step.   Equation 
2.99  defines  parameter  <    as  the  inverse  of  parameter  k.   The 
variation  of  k  with  flow  is  defined  by  equation  4.35. 

Equation  2.84  is  still  valid  for  nonlinear  reservoirs 
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£!■«  (i-p)«i  +p-i  -  (i-p)-o  -p-o 

At  i       2  1       2 

Rewriting  AS  in  terms  of  storage  at  the  beginning  and  end  of  a  time  step 
for  the  itn  reservoir  yields 

S   .  .  -  S 


t±4_   t'1-di»HttH+p.ottltl.1- 


(4.37) 


a-p)-°t,i-p*Vi,i 


Substituting  for  S. }i  with  equation  4.36  produces 


t+l,i   t+1,1    t,i  t,i  = 
At 

(1^)-°t,i-i  +  p-Vi,i-i-  (1^)*°t,i-p*Vi,i 


(4.38) 


Solving  for  Ot+isi 

[<          -   (l-p)-At) 
o  (i-P)-^       ,Q  +     l   t,i        ^_.0     _   + 

6+1,1        C^,,  +  P-Atj      t'1"1        (<t+1,i+P'At) 

(4.39) 

7  ~T   t+i,i-i 

(ict+11  +  P'At) 


The  value  <t+1  i  in  equation  4.39  depends  upon  Ot+1,1-   To  avoid  the 
need  for  an  iterative  solution  let 

k   ,  .  -  <  4  (4.40) 

t+l,i    t,i 

Using  the  approximation  of  equation  4.40  and  defining 

E    =  k   ,  +  P'At  (4.41) 

t,i    t,i 
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and 


(*tfl  "  (l-P>At] 


1  . t>x  , (4.42) 

t'1         5t,i 


equation  4.39  becomes 

(l-p).At  ,  +^0   ,,,       (4.43) 

t+i.i       i„  .       t,i-i      t,i   t,i    r~7  t+ifi-i 

t  ,  1  C  ,  1 

For  the  first  reservoir  in  the  cascade,  0.  j_j  is  the  specified  system 

inflow,  so  equation  4.43  reduces  to 

Af((l-p)«I   +  P'It+1) 

0   ,  ,  -*   •()  ,  +  ' ElL-  (4.44) 

t+1,1    t,l   t,l  4t  t 

By  defining 

Vi  =  (^>it+P'Vi  (4-45) 

the  governing  equation  for  the  first  reservoir  in  the  cascade  becomes 

0      =  *   -0   ,  +  ,At  'u  ,  (4.46) 

t+i,i     t,i  t,i    ftt  t+i 
t ,  1 

The  governing  equations  for  subsequent  reservoirs  can  be  found  by 

recursively  substituting  for  0t+1>1_!  in  equation  4.43. 

Define 

T     =  *L_  (4.47) 

t«1    ^7 


and 


fl   .  =  (1-p)  +  p»t   .  (4.48) 

t ,  i  t  ,i 


The  value  of  the  element  for  the  jtn  row  and  mtn  column  of  the  j^  t 
matrix  in  equation  1.1  is 
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zero  for  m  >  j  (4.49) 

¥  for  m  =  j  (4.50) 

t  ,tn 

fl   -T  for  m  =  1-1  (4.51) 

t,m   t,j 

a        .p^"™"1.   n   T  t  for  in  <  j-1  (4.52) 

**"  Wl  t'1 

where  II  is  the  product  operator. 


The 


ith 


jcri  row  of  the  G^   fc  matrix  i 

j 


J"1,  n  T   .  (4.53) 


1=1 


t,i 


Equations  4.49  through  4.53  complete  the  information  needed  to 
construct  a  nonlinear  state  space  routing  model.   The  steps  that  have 
led  to  this  point  can  he  summarized  as 

1.  development  of  a  linear  state  space  routing  model,  having 
constant  parameters  a,  n  and  k,  with  an  impulse  response 
identical  to  the  linearized  Saint-Venant  equations, 

2.  determination  of  the  variations  with  flow  level  of  the  mean 
and  variance  of  the  linear  state  space  model  impulse 
response  function, 

3.  relation  of  the  variations  with  flow  of  the  impulse  response 
function's  moments  to  variations  of  the  parameters  a,  n  and 
k,  and 

4.  derivation  of  the  state  transition  and  input  coefficient 
matrices  as  functions  of  flow  level. 

The  theory  is  now  well  defined:   let  the  model  parameters  vary 
with  flow  and  simulate  using  equation  1.1.   There  are,  however,  many 
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ways  to  implement  the  theory,  and  many  possible  algorithms  for 
simulation.   The  general  approach  in  developing  the  nonlinear  state 
space  model  is  to  retain  as  much  of  the  theory  as  possible  while  keeping 
the  model  algorithm  simple.   The  final  algorithm,  presented  here, 
incorporates  all  the  theory  developed  in  this  dissertation.   The 
preliminary  algorithms  are  described  in  the  following  section.  As  more 
of  the  theory  was  introduced  into  the  approaches  to  state  space 
modelling  of  unsteady  flow,  the  results  compared  more  favorably  with  the 
solution  of  the  full  Saint-Venant  equations.   The  final  algorithm 
emulates  well  the  solution  of  the  complete  equations  of  unsteady  flow 
over  a  wide  range  of  channel  and  wave  conditions. 

The  algorithm  for  simulating  with  the  nonlinear  state  space  model 


is: 


1.  Compute  the  average  of  the  outflows  from  each  of  the  n 
reservoirs  at  time  =  t,  and  the  current  system  inflow 

Q    -  (u   ,  +  I  0   J/(n+l)  (4.54) 

t+l    t+i   /\  t,r 

where  I    is  the  summation  operator. 

2.  Use  7)         with  equation  4.33  and  nQ  at  q   (the  number  of 
reservoirs  at  the  reference  flow  level)  to  find  the  number 
of  reservoirs  needed  (nt+1>.  and  round  this  value  to  the 
nearest  integer. 

3.  If  nt+j  ^  nt,  interpolate  the  reservoir  outflows  to  compute 
a  new  Xj-  state  vector  of  dimension  nt+^. 

4.  Compute  the  average  flow  through  each  reservoir; 
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for  the  first  reservoir 

q    =  (u    +0  J/2  (4.55) 

t,l   k  t+1    t,\ 

and  for  subsequent  reservoirs 

°t,i  -  (°,.i  +  °t,i-i)/2  <4"56) 

5.  Compute  values  for  at  .  and  kt  i   for  the  flow  levels 
determined  in  step  4,  using  equations  4.31  and  4.35,  and  the 
values  for  aQ  and  kQ  at  Q  Q . 

6.  Compute  parameter  a  as  the  average  of  the  delays  for  the 
flow  through  each  of  the  nt+j  reservoirs,  as 


I  a  J/n 
i=l  C,:L 


(4.57) 


7.  Compute  the  elements  of  _F^  t  and  Gj(  t  with  equations  4.49 
through  4.53. 

8.  Compute  ^+^   with  equation  1.1. 

9.  Lag  the  output  by  ¥   ,  and 

10.   increment  the  time  and  return  to  step  1. 

This  algorithm  allows  all  parameters  of  the  nonlinear  state  space 
model  to  vary  with  flow  as  specified  by  equations  4.31,  4.33  and  4.35. 
The  thesis  of  this  study  is  that  these  relationships  can  emulate  the 
behavior  of  the  full  Saint-Venant  equations.   Results  verifying  this 
thesis  are  presented  in  Chapter  V. 

The  Path  to  the  Final  Model 


The  algorithm  for  the  final  nonlinear  state  space  routing  model 
performs  the  necessary  sequence  of  steps  to  assure  that  the  variations 
with  flow  of  the  moments  of  the  impulse  response  function  for  the 
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linearized  Saint-Venant  equations  obey  the  previously  derived 
theoretical  relationships.  The  approach  taken  to  reach  that  final 
algorithm  was  to  begin  with  a  simple  algorithm  and  add  complexity  only 
if  it  improves  comparisons  with  the  full  unsteady  flow  model. 

The  major  factor  complicating  the  implementation  of  the  nonlinear 
state  space  model  is  the  variation  with  flow  of  the  number  of 
reservoirs,  n.   There  are  n  elements  in  the  state  vector,  X_.   Letting  n 
vary  with  flow  implies  that  the  dimensionality  of  X_  changes  over  time. 
Initial  algorithms  were  designed  to  keep  n  fixed.   That  simplified  the 
algorithm,  but  immediately  introduced  theoretical  problems. 

The  lag  through  a  cascade  of  linear  reservoirs  is 

(4.58) 

Equations  4.31,  4.33  and  4.35,  which  define  the  variations  with  flow  of 
a,  n  and  k,  can  be  substituted  into  equation  4.58  to  produce 


((r-l)/m 

h   ~  Q'~      +  „(m-2+r)/ 


(l-m)/m  +  _0_____  (4.59) 


_  Q(l-m)/m  +  Q(l-m)/m  (4>6Q) 

Both  terms  on  the  right  hand  side  of  equation  4.60  vary  with  flow  to  the 
(l-m)/m  power.   If  equation  4.33  is  simplified  as 

n  =  n  (^-61) 

0 

where   n   =  a  constant,  then,  in  order  to  maintain  the  variation  with 
flow  of  the  lag  indicated  by  equation  4.60,  the  variation  of  k  must 
become 
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k  ~  Q 


(m-l)/m 


(4.62) 


Equations  4.62  and  4.35  are  not  consistent.   Holding  the  number  of 
reservoirs  constant  warps  the  variation  with  flow  of  k.   Because  k  is 
also  the  inverse  of  K  in  equation  4.23,  the  storage  computed  for  each 
reservoir  is  no  longer  that  predicted  by  the  theory. 

If  the  pure  delay,  a,  is  held  constant,  equation  4.60  becomes 


t£  =  Q  +  Q 


0   ^(l-m)/m 


(4.63) 


which  is  inconsistent  with  previous  results.   The  variation  of  parameter 
k  cannot  obey  both  equations  4.35  and  4.63. 

In  spite  of  these  theoretical  shortcomings,  several  simpler 
algorithms  were  developed  for  a  nonlinear  state  space  routing  model.   In 
each  case,  the  simplification  led  to  a  theoretical  inconsistency.   The 
initial  algorithms,  with  n  =  nQ ,  are: 

A)   1.   Fix  a  =  a  at  a  reference  flow  level. 

2.  Compute  the  lag  for  each  reservoir  with  equation  4.7. 

3.  Solve 


t.i 


(4.64) 


(t, 


't.i 


4.  Compute  j^  t  and  _G_x  t  using  the  kt  j_'s  and  equations 
4.49  through  4.53. 

5.  Predict  _Xt+l  with  equation  1.1,  and 

6.  Delay  the  output  by  a$  • 
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B)   1.  Fix  k  =  1c,  at  a  reference  flow  level  (Compute  J_  and  G_ 
for  k  =  k  ). 

2.  Compute  the  lag  for  the  average  flow  through  the  cascade 
with  equation  4.7. 

3.  Solve 

n 

a  =  t   -   0  (4.65) 

t    JL    k 

t    o 

4.  Predict  Xr+1  with  equation  1.1,  and 

5.  Delay  the  output  by  at . 


C)   1.   Compute  kt  *    f°r  each  reservoir  with  equation  4.7. 

2.  Compute  the  lag,  to    for  the  average  flow  through  the 

*t 

cascade  with  equation  4.7. 

3.  Solve 

n 

«t-H  -1°  r-  (4*66) 

C    *t   i=l   t,i 

4.  Compute _Fy  t  and  ^  t  using  the  kt  ^'s  and  equations 
4.49  through  4.53 

5.  Predict  2^.+1  with  equation  1.1,  and 

6.  Delay  the  output  by  at • 

Theoretical  inconsistencies  caused  those  preliminary  algorithms  to 
emulate  poorly  the  behavior  of  the  primary  model.   Parameter  variations 
with  flow  were  warped  by  the  improper  governing  equations. 

The  next  step  on  the  path  to  the  final  algorithm  was  to  let  the 
number  of  reservoirs  vary  with  flow.   That  allowed  theoretical 
consistency  to  be  maintained,  but  incresed  the  complexity  of  the 
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algorithm.   At  this  point  the  number  of  state  variables  could  change 
over  time;  therefore,  some  mechanism  had  to  be  developed  to  produce  a 
new  state  vector  Xj-  with  nt+i  states  from  an  existing  Xj-  with  nt 
states.   All  the  algorithms  studied  thus  far  had  assumed  that  the  states 
(the  elements  of  vector  X)  were  reservoir  storages.   Reservoir  storage 
and  outflow  are  interchangable,  given  kt  if  with  equation  4.36. 
However,  when  the  number  of  states  can  vary,  storage  is  an  inconvenient 
choice  for  the  state  variable.   The  approach  taken  to  change  the  number 
of  states  when  the  states  are  storage  is  to  maintain  continuity.   That 
is,  to  keep  the  same  total  volume  of  water  in  the  nt+^  reservoirs  as  is 
originally  in  the  nt  reservoirs. 

Letting  n  vary  when  the  states  are  storage  created  problems 
because  the  value  of  kt>i  changes  smoothly,  while  nt+1  is  an  abrupt 
change  from  ru.   The  simulated  outflow  for  the  channel  (i.e.,  from  the 
nt+,  reservoir),  computed  with  equation  4.36,  has  discontinuities  when  n 
changes.   Those  jumps  in  the  flow  become  more  exaggerated  as  n  becomes 
small,  because  the  percentage  change  between  nt  and  nt+^  increases. 

To  circumvent  those  problems,  the  algorithm  was  rewritten  with 
reservoir  outflow  as  the  states.  When  n  changes,  the  flows  are  linearly 
interpolated  to  produce  the  new  states,  and  no  abrupt  changes  occur 
because  the  downstream  outflow  is  fixed  during  the  interpolation. 

The  algorithm  is  now  essentially  in  final  form.   Interpretation  of 
the  flow  levels  at  which  to  compute  the  parameters,  n,  a  and  k,  is  the 
last  step.   Originally,  the  parameters  of  the  nonlinear  state  space 
model  varied  with  the  reservoir  outflow.   The  new  number  of  states, 
nt+l,  was  found  with  equation  4.33  using 


v.  ■  J,  °t,i 


158 


(4.67) 


Th 


e  at  i's  and  kt  j's  were  found  with  equations  4.31  and  4.35  using 

0    =0  (4.68) 

t,i    t,i 

Equations  4.67  and  4.68  look  only  at  the  reservoir  outflows  to  solve  for 
a,  n  and  k.   The  equations  contained  in  the  final  algorithm  (equations 
4.54  through  4.56)  use  an  average  of  the  inflow  and  outflow  to  better 
represent  the  conditions  in  the  reservoir. 

These  steps  lead,  finally,  to  the  nonlinear  state  space  routing 
algorithm  presented.   All  the  parameters  vary  with  flow  level.   In  this 
way  the  properties  of  the  Saint-Venant  equations'  impulse  response 
function  are  matched  by  the  nonlinear  state  space  model. 

Summary 

A  nonlinear  state  space  model  that  can  emulate  the  behavior  of  the 
equations  of  unsteady  flow  has  been  developed.   Results  obtained  by 
simulating  with  both  the  nonlinear  state  space  model  and  the  DWOPER 
model's  numerical  solution  of  the  Saint-Venant  equations  are  presented 
in  Chapter  V. 

In  general,  however,  any  surface  runoff  or  channel  routing  model 
can  be  chosen  as  the  primary  model.   The  strategy  employed  to  emulate 
any  primary  model  is  shown  in  Figure  4.12.   The  steps  needed  to  fulfill 
this  strategy  are  presented  in  Chapter  VI. 
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NONLINEAR  STATE  SPACE  MODEL 
AS  SURROGATE  FOR  PRIMARY  MODEL 


Figure  4.12    General  Strategy  to  Emulate 
Any  Primary  Model 


CHAPTER  V 
LIMITATIONS  OF  THE  NONLINEAR  STATE  SPACE  APPROACH 


Introduction 


The  nonlinear  state  space  model  developed  in  Chapter  IV  models  the 
stationary  nonlinear  behavior  of  the  full  Saint-Venant  equations  with  a 
nonstationary  linear  structure.   Figure  5.1  shows  the  process  used  to 
verify  that  the  nonlinear  state  space  model  can  emulate  the  equations  of 
unsteady  flow.   Given  a  prismatic  channel,  the  value  of  m  can  be  found 
with  equation  2.28.   For  prismatic  channel  shapes  defined  by  equation 
2.23,  the  value  of  m  will  be  the  same  for  all  flows.  However,  in 
general,  m  can  vary  with  flow  if  the  channel  shape  if  irregular.   This 
possibility  is  explored  in  Chapter  VI. 

For  results  presented  in  this  chapter,  m  is  constant  for  all  flows 
and  can  be  found  directly  from  equation  2.28.   Given  m,  the  variations 
with  flow  of  parameters  a,  n  and  k  are  specified  by  equations  4.31,  4.33 
and  4.35.   Reference  points  to  fix  the  parameter  values  for  a  particular 
channel  shape,  flow  level  and  distance  from  the  channel  inlet  are  given 
in  Table  3.4.   Recall  that  these  parameter  values  were  computed  from  the 
moments  of  the  inflow  and  outflow  hydrographs  by  first  using  equations 
2.53  and  2.54  to  find  the  mean  and  variance  of  the  system  response 
function,  and  then  applying  equations  2.78  through  2.80.   Results  from 
the  state  space  and  OWOPER  models  for  large  inflow  transients  are 
compared  in  the  next  section. 
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Figure  5.1    Verification  Process  for  the 
Nonlinear  State  Space  Model 
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Verification  of  the  Nonlinear  State  Space  Model 
Introduction 

The  nonlinear  state  space  and  DWOPER  models  are  compared  for 
various  inflow  hydrographs  in  rectangular  and  triangular  channels 
described  in  Chapter  III.   The  remaining  figures  in  this  chapter  are 
organized  as  shown  in  Table  5.1.   Each  of  the  16  figures  used  to  verify 
the  nonlinear  state  space  model  corresponds  to  a  figure  in  Chapter  III, 
comparing  the  linear  state  space  and  DWOPER  models.   The  pairs  of 
corresponding  figures  are  listed  in  Table  5.2. 

Results  for  several  single-peaked  inflow  hydrographs  and  a 
double-peaked  hydrograph  are  discussed  in  the  following  sections.   In 
addition,  the  issues  of  continuity,  the  number  of  reservoirs  required  by 
the  state  space  model,  and  the  CPU  times  to  simulate  these  results  are 
presented. 

Single-Peaked  Inflow  Hydrographs 

A  very  rapid  transient  (<?in  =  3  hours)  routed  in  a  rectangular 
channel  for  100  miles  is  shown  in  Figure  5.2.   The  nonlinear  state  space 
model  does  not  properly  lag  the  rising  limb  of  the  outflow  hydrograph. 
This  problem  appears  to  be  almost  the  same  as  that  encountered  for  the 
linear  state  space  model  in  Figure  3.28.   The  causes  of  these  problems 
are  related,  but  the  differences  are  important  to  understanding  the 
nonlinear  state  space  model.   In  Chapter  III,  the  parameters  of  the 
state  space  model  were  fixed  at  the  values  for  the  reference  flow,  00 • 
In  Figure  3.28,  QQ  equaled  100000  cfs.  A  very  peaked  inflow  hydrograph, 
as  when  a±  =3  hours,  disperses  rapidly  as  it  moves  downstream.   The 
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Table  5.1   Organization  of  Figures  5.2  through  5.17, 
Comparing  the  Nonlinear  State  Space  Model 
with  the  Full  Saint-Venant  Equations 
for  Large  Inflow  Transients 


Rectangular  Channel 


L  =  100  miles 


L  =  400  miles 


oin  =  3  hours 
a.  =  6  hours 
o^n  =  12  hours 
a.      =  24  hours 


Figure  5.2 
Figure  5.4 
Figure  5.6 
Figure  5.8 


Figure  5.3 
Figure  5.5 
Figure  5.7 
Figure  5.9 


Triangular  Channel 


°in 


=   6  hours 


=  12  hours 


L  =  100  miles 

Figure  5.10 
Figure  5.12 


L  =  400  miles 

Figure  5.11 
Figure  5.13 


Double-Peaked  Hydrograph 


00  =  100000  cfs 
QQ  =  10000  cfs 
QQ  =  200000  cfs 


L  =  100  miles 
Figure  5.14 


L  =  400  miles 

Figure  5.15 
Figure  5.16 
Figure  5.17 
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Table  5.2   Pairs  of  Figures  with  Identical  Input  Data 

Comparing  the  Linear  and  Nonlinear  State  Space 
Models  with  the  Full  Saint-Venant  Equations 
for  Large  Inflow  Transients 


Figure  Comparing  Linear        Figure  Comparing  Nonlinear 
State  Space  Model  and  Full       State  Space  Model  and  Full 
Saint-Venant  Equations  Saint-Venant  Equations 


3.28  5.2 

3.29  5.3 

3.30  5.4 

3.31  5.5 

3.32  5.6 

3.33  5.7 

3.34  5.8 

3.35  5.9 

3.36  5.10 

3.37  5.11 

3.38  5.12 

3.39  5.13 

3.40  5.14 

3.41  5.15 

3.42  5.16 

3.43  5.17 
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peak  flow  a  short  distance  from  the  channel  inlet  is  significantly  lower 
that  the  inflow  peak  of  210000  cfs.   This  caused  problems  for  the  linear 
state  space  model  when  its  parameters  were  selected  for  QQ  =  100000  cfs, 
because  the  flow  in  the  channel  was  below  this  level  the  vast  majority 
of  the  time.   As  we  have  seen  in  Chapter  IV,  channel  lag  decreases  as 
the  flow  level  increases.   This  means  that,  since  the  actual  flow  was 
almost  always  less  than  the  reference  flow  for  the  linear  state  space 
model,  it  did  not  delay  the  hydrograph  sufficiently.   Therefore,  the 
entire  outflow  hydrograph  simulated  by  the  linear  state  space  model  in 
Figure  3.28  was  earlier  than  the  DWOPER  model  results. 

The  nonlinear  state  space  model  does  not  adequately  lag  the  rising 
limb  of  the  inflow  hydrograph  for  a  related  but  different  reason.   The 
parameters  of  the  nonlinear  state  space  model  are  functions  of  the 
average  flow  level  in  the  channel.   The  average  flow  in  the  first 
reservoir  is  defined  by  equation  4.55  as 

When  an  inflow  hydrograph  rises  very  steeply,  this  equation  gives  false 
information  about  the  actual  average  flow  for  that  reservoir.   The 
assumption  implicit  in  equation  4.55  is  that  the  average  flow  is  equal 
to  the  average  of  the  reservior  inflow  and  outflow.   If  the  inflow 
disperses  very  rapidly,  equation  4.55  provides  a  poor  representation  of 
the  actual  average  flow.   The  average  flow  predicted  by  equation  4.55 
will  always  be  high  for  the  rising  limb  of  a  hydrograph.   For  a  very 
rapidly  rising  hydrograph  it  will  be  significantly  high.   For  a  short 
channel  length  the  problem  is  compounded  because  the  total  number  of 
reservoirs  is  small  and,  therefore,  the  first  reservior  has  a  larger 
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impact  on  the  outflow  hydrograph.   The  parameters  of  the  first  reservoir 
of  the  nonlinear  state  space  model  are  computed  at  a  reference  flow 
level  higher  than  actually  in  the  channel.   Therefore,  as  with  the 
linear  state  space  model,  the  outflow  is  earlier  than  the  DWOPER  model 
results . 

There  is,  however,  a  substantial  improvement  in  the  recession,  or 
falling  limb,  of  the  state  space  model  outflow  hydrograph  from  Figure 
3.28  to  Figure  5.2.   Once  the  rapidly  varying  portion  of  the  inflow 
hydrograph  has  entered  the  channel,  the  equations  for  choosing  reference 
flow  levels  are  much  more  representative  of  the  actual  channel 
conditions.   Since  the  lag  and  dispersive  properties  of  the  nonlinear 
state  space  model  can  vary  by  changing  parameters  values,  the  model  is 
able  to  match  the  recession  portion  of  the  DWOPER  model  outflow 
hydrograph. 

At  a  distance  of  400  miles  from  the  channel  inlet,  shown  in  Figure 
5.3,  results  from  the  nonlinear  state  space  model  are  greatly  improved 
over  those  of  the  linear  state  space  model  in  Figure  3.29.   The  outflow 
hydrograph  from  the  linear  model  arrives  much  too  early  because  its 
parameters  are  fixed  at  0   =  100000  cfs.   The  nonlinear  state  space 
model  is  much  closer  to  the  lag  of  the  DWOPER  model  because  its 
parameters  vary  with  flow  level  in  the  channel.   The  rising  limb  of  the 
nonlinear  state  space  model  outflow  hydrograph  is  too  early  in  Figure 
5.3  for  the  same  reasons  given  for  Figure  5.2.   The  situation  is 
ameliorated  at  400  miles  from  the  channel  inlet,  however,  because  the 
total  number  of  reservoirs  is  increased  and  the  first  reservoir  has  a 
proportionally  smaller  influence  on  the  model  outflow. 
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Figure  5.4  presents  results  for  a  less  abrupt  inflow  hydrograph 
(ain  =  6  hours)  at  L  =  100  miles.  Again,  there  are  some  discrepancies 
between  the  nonlinear  state  space  and  DWOPER  models  for  the  rising  limb 
of  the  hydrograph  which  are  caused  by  an  improper  representation  of  the 
average  flow  in  each  reservoir  and,  therefore,  an  incorrect  choice  of 
state  space  model  parameter  values.   After  the  abrupt  portion  of  the 
inflow  hydrograph  has  entered  the  channel,  the  nonlinear  state  space 
model  agrees  very  closely  with  the  recession  portion  of  the  DWOPER  model 
outflow  hydrograph.   This  is  a  vast  improvement  over  the  linear  model  in 
Figure  3.30,  which  displays  the  almost  symmetrical  shape  typical  of  a 
linear  model. 

A  a.      =  6  hours  inflow  hydrograph  routed  400  miles  in  a 
rectangular  channel  is  shown  in  Figure  5.5.   The  timing  of  the  peak 
predicted  by  the  state  space  model  is  very  close  to  that  from  the  DWOPER 
model.   Although  the  rising  limb  arrives  at  the  channel  outlet  too  soon, 
the  shape  of  the  recession  is  a  significant  improvement  over  the  linear 
model  in  Figure  3.31. 

A  comparison  of  Figures  3.32  and  5.6  provides  a  classic  example  of 
a  linear  versus  a  nonlinear  model  response.   The  nonlinear  state  space 
model  for  oin  =  12  hours  at  100  miles  in  Figure  5.6  produces  an  outflow 
hydrograph  with  sharply  rising  and  slowly  falling  limbs,  and  matches  the 
DWOPER  model  response  very  well.   Outflow  from  the  linear  state  space 
model  in  Figure  3.32  rises  too  slowly  and  falls  to  rapidly.   It  is 
doomed  by  its  linear  structure  to  fail  when  trying  to  emulate  a 
nonlinear  model. 

At  L  =  400  miles,  in  Figures  3.33  and  5.7,  the  linear  versus 
nonlinear  effect  is  even  more  pronounced.   Although  the  outflow  from  the 
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nonlinear  state  space  model  begins  to  rise  too  soon,  the  steepness  of 
its  ascent  is,  in  general,  equal  to  that  of  the  DWOPER  model.   This  is 
in  sharp  contrast  with  the  Gaussian  shape  of  the  linear  model  in  Figure 
3.33. 

The  linear  state  space  model  performed  adequately  for  the  more 
slowly  varying  inflow  hydrograph  (oin  =  24  hours)  at  100  and  400  miles 
from  the  channel  inlet  as  was  seen  in  Figures  3.34  and  3.35.   As  shown 
in  Figure  5.8,  the  major  improvements  introduced  by  the  nonlinear  state 
space  model  at  100  miles  are  at  the  beginning  of  the  rising  and  end  of 
the  recession  limbs  of  the  outflow  hydrograph.   The  DWOPER  model  is 
matched  very  well  by  the  nonlinear  state  space  model.   At  400  miles  from 
the  channel  inlet,  shown  in  Figure  5.9,  the  nonlinear  state  space  model 
has  produced  an  outflow  hydrograph  closer  to  that  from  the  DWOPER  model 
than  was  possible  with  the  linear  model. 

As  mentioned  in  Chapter  III,  the  incremental  improvement  from 
using  a  linear  rather  than  a  nonlinear  model  would  be  expected  to  be 
greater  for  a  rectangular  channel  than  for  a  triangular  channel.   The 
behavior  of  a  triangular  channel  is  more  linear  than  for  a  rectangular 
channel.   However,  even  for  a  triangular  channel,  the  nonlinear  state 
space  model  shows  improvement  over  the  linear  state  space  model. 

Comparison  of  Figures  3.36  and  5.10  for  a^n  =  f,   hours  in  a 
triangular  channel  at  L  =  100  miles  shows  that  the  nonlinear  state  space 
model  approximates  the  recession  portion  of  the  DWOPER  model  outflow 
hydrograph  better  than  the  linear  model.   At  400  miles  from  the  channel 
inlet,  as  shown  in  Figure  5.11,  the  nonlinear  state  space  model  outflow 
closely  resembles  the  DWOPER  model  results  except  for  an  early  rising 
limb.   The  timing  of  the  peak  flow  is  much  improved  over  the  linear 
model  outflow  seen  in  Figure  3.37. 
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For  04   =  12  hours  in  a  triangular  channel,  results  from  the 
nonlinear  state  space  model  compare  well  with  the  DWOPER  model,  as  shown 
in  Figures  5.12  and  5.13.   Although  the  timing  of  the  linear  state  space 
model,  seen  in  Figures  3.38  and  3.39,  was  good,  the  nonlinearities  in 
the  rising  and  falling  limbs  are  modelled  significantly  better  by  the 
nonlinear  state  space  model. 

Typical  hydrograph  durations  for  natural  channels  vary  widely 
(Chow  1964,  p.  25-5).   An  empirical  relationshop  for  the  time  base  of  a 
hydrograph  was  developed  by  Snyder  (1938)  as 
t 


T  =  3  + 


(5.1) 


where   t  is  the  lag  time  to  the  peak  of  a  unit  hydrograph  in  hours,  and 

T  has  units  of  days. 
The  unit  hydrograph  peak  lag  time  can  be  expressed  as 
L«L 


t   =  C 


(5.2) 


P    C   /I 

where   Ct  ±s   an  empirical  coefficient  (usually  equal  to  0.35  for  a 
valley  drainage  area, 
L   is  the  river  mileage  through  the  basin, 
L   is  the  river  mileage  from  the  basin  outlet  to  the  basin 

center  of  mass, 
s    is  the  basin  slope  in  feet  per  mile,  and 

n   is  an  empirical  coefficient  equal  to  0.38  (Linsley,  Kohler 
and  Paulhus  1958,  p.  207). 
An  approximation  of  the  hydrograph  duration  can  be  found  with  equations 
5.1  and  5.2.   Hydrographs  with  o^n  >  12  hours  would  be  expected  for  most 
basins.   For  inflow  hydrographs  of  this  duration  equation  4.55 
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adequately  captures  the  rising  limb,  and  the  nonlinear  state  space  model 
emulates  well  the  primary  DWOPER  model. 

A  Double-Peaked  Inflow  Hydrograph 

At  100  miles  from  the  channel  inlet,  the  nonlinear  state  space 
model  approximates  well  the  DWOPER  model  for  the  double-peaked  inflow 
hydrograph  shown  in  Figure  5.14.   As  with  results  presented  above,  the 
major  improvement  of  the  nonlinear  model  over  the  linear  state  space 
model  is  in  the  shape  of  the  rising  limb  and  in  the  recession  of  the 
outflow  hydrograph. 

In  Figure  5.15,  a  double-peaked  hydrograph  is  routed  400  miles 
with  the  nonlinear  state  space  and  DWOPER  models.   Results  from  the  two 
models  resemble  each  other  very  well  for  the  entire  outflow  hydrograph 
duration.   Figures  5.16  and  5.17  look  identical  to  Figure  5.15. 
Nonlinear  state  space  model  parameters  were  chosen  at  Q  =  10000  and 
200000  cfs  in  Figures  5.16  and  5.17,  respectively.   These  figures  should 
be  compared  with  Figures  3.42  and  3.42  for  the  linear  state  space 
model.   The  reference  flow  level  at  which  parameters  are  selected  is 
critical  to  the  linear  state  space  model  behavior.   The  variations  with 
flow  level  expressed  in  equations  4.31,  4.33  and  4.35  make  the  nonlinear 
state  space  model  insensitive  to  the  flow  level  at  which  its  reference 
parameter  values  are  chosen. 

Continuity 

The  linear  state  space  model  was  guaranteed  to  preserve  continuity 
(i.e.,  to  make  certain  that  the  water  leaving  the  channel  equaled  the 
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water  entering).  With  constant  parameters,  each  reservoir  in  the  linear 
cascade  let  all  the  water  that  entered  pass  through  with  no  increase  or 
decrease  in  quantity.   However,  there  is  no  such  guarantee  with  the 
nonlinear  state  space  model.   In  fact,  it  is  almost  certain  that  the 
volume  of  an  inflow  transient  will  be  modified  in  some  way  because  of 
the  variations  in  parameter  values. 

As  with  any  discrete-in-time  model,  the  nonlinear  state  space 
model  produces  an  outflow  value  each  time  step.   The  pure  delay,  a, 
changes  with  the  flow  level  in  the  channel.   The  time  for  which  the 
outflow  applies  is  equal  to  the  time  steps  taken  since  the  simulation 
began  plus  an  adjustment  for  the  pure  delay.   This  means  that  the  time 
period  for  which  an  outflow  value  applies  may  not  be  one  time  step  after 
the  previous  value  simulated.  When  flow  in  the  channel  is  increasing, 
the  value  of  a  will  decrease.   This  allows  the  nonlinear  state  space 
model  to  approximate  the  steep  rising  limb  of  the  outflow  hydrograph  of 
the  DWOPER  model.   Similarly,  the  slow  decay  of  the  hydrograph' s 
recession  is  modelled  by  increasing  the  pure  delay  as  flow  in  the 
channel  decreases.   However,  this  means  that  the  outflow  from  the  last 
reservoir  may  apply  for  less  or  more  than  a  single  time  step.   During 
the  course  of  an  inflow  transient  and  return  to  baseflow  level,  the 
errors  in  continuity  should  compensate. 

Results  presented  in  Table  5.3  indicate  the  degree  to  which  this 
is  true.   The  steeply  rising  inflow  hydrographs  (oin  =  3  hours)  cause 
the  nonlinear  state  space  model  to  diverge  from  continuity  by  7  to  9 
percent.   For  all  other  inflows,  the  nonlinear  state  space  model  is 
within  plus  or  minus  3.3  percent  of  maintaining  continuity. 


Table  5.3   Continuity  Check  for  the  Nonlinear  State 
Space  and  DWOPER  Models  for  Results  of 
Figures  5.2  through  5.17 


Figure 

Ratio  of  Outflow 

to  Inflow  for  the 

Nonlinear  State 

State  Model 

Ratio  of  Outflow 

to  Inflow  for  the 

DWOPER  Model 

5.2 

1.070 

0.999 

5.3 

1.093 

0.985 

5.4 

0.999 

1.000 

5.5 

1.033 

0.985 

5.6 

0.995 

1.000 

5.7 

0.997 

1.000 

5.8 

0.998 

1.000 

5.9 

0.995 

1.000 

5.10 

1.026 

1.000 

5.11 

1.030 

0.999 

5.12 

0.998 

1.000 

5.13 

1.003 

0.999 

5.14 

0.995 

1.000 

5.15 

0.988 

1.000 

5.16 

0.991 

1.000 

5.17 

0.990 

1.000 
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The  Number  of  Reservoirs 

Conceptually,  the  DWOPER  model  Is  similar  to  a  cascade  of 
nonlinear  reservoirs.   Each  of  the  distance  steps  along  the  channel  at 
which  the  Saint-Venant  equations  are  solved  can  be  thought  of  as  a 
nonlinear  reservoir  governed  not  by  equation  4.23,  but  by  the  equations 
of  unsteady  flow.   The  number  of  distance  steps  used  to  simulate  the 
results  presented  above  ranged  from  10  for  a^n  =  24  hours  to  160  for 
ain  =  ^  hours.   The  number  of  reservoirs  in  the  nonlinear  state  space 
model  is  a  function  of  flow  level  in  the  channel.   It  must  be  an  integer 
each  time  step,  but  varies  over  time  as  the  inflow  transient  rises  and 
ebbs.   The  average  numbers  of  reservoirs  in  the  nonlinear  state  space 
model  cascade  are  presented  in  Table  5.4.   The  ratio  of  DWOPER  model 
distance  steps  to  nonlinear  state  space  model  reservoirs  varies  from  1:1 
to  13:1.   We  would  expect  the  average  number  of  reservoirs  to  decrease 
as  the  o-!      value  increases  because  the  flow  level  in  the  channel  is 
higher  for  a  longer  period  of  time  and,  therefore,  the  number  of 
reservoirs  decreases.   This  is  not  seen  in  Table  5.4  because  the  total 
simulation  period  was  increased  as  a^n  increased  to  capture  the  entire 
outflow  transient. 

CPU  Time 


Although  differences  in  the  nonlinear  state  space  and  DWOPER  model 
results  shown  in  Figures  5.2  through  5.17  are  generally  small, 
differences  in  the  CPU  time  required  are  sometimes  large.   The  CPU  times 
required  by  all  the  nonlinear  state  space  and  DWOPER  model  sumulations 
are  presented  in  Table  5.5.   The  DWOPER  model  requires  approximately  2 
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Table    5.4      The   Average   Number  of   Reservoirs   Used 
by   the  Nonlinear  State   Space  Model 
in  Figures    5.2  through    5.17 


Rectangular  Channel 


a,  =  3  hours 
°in  =  6  hours 
ain  =  ^  hours 
ain  =   24  hours 


L   =    100  miles 

3.533 
3.078 
3.085 
3.201 


L  =  400  miles 

12.212 
10.375 
11.050 
11.894 


Triangular  Channel 


L  =  100  miles 


L  =  400  miles 


'in 


2.712 

2.413 


9.571 
8.534 


Double-Peaked  Hydrograph 


0Q  =  100000  cfs 
Q0  =  10000  cfs 
Q„  =  200000  cfs 


L  =  100  miles 
3.057 


L  =  400  miles 

11.227 
11.140 
11.243 
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Table  5.5  CPU  Times  for  the  Nonlinear  State  Space  and 
DWOPER  Models  on  a  PRIME  7  50  Computer  for 
Results  of  Figures  5.2  through  5.17 


Rectangular  Channel 


Nonlinear  State 

Figure 

,   °in  % 
(hours) 

3 

L 
(miles) 

100 

% 
(1000  cfs) 

100 

Space 
CPU 

Model 
Time 

DWOPER  Model 
CPU  Time 

5.2 

1.31 

sec 

— 

5.3 

3 

400 

100 

5.15 

sec 

409.06  sec 

5.4 

6 

100 

100 

0.27 

sec 

— 

5.5 

6 

400 

100 

0.74 

sec 

175.87  sec 

5.6 

12 

100 

100 

0.23 

sec 

— 

5.7 

12 

400 

100 

0.63 

sec 

63.80  sec 

5.8 

24 

100 

100 

0.35 

sec 

— 

5.9 

24 

400 

100 

0.95 

sec 

23.58  sec 

Triangular  Channel 


Figure 

,  ain  . 
(hours) 

6 

L 
(miles) 

100 

% 
(1000  cfs) 

100 

Nc 

nlinear  State 
Space  Model 
CPU  Time 

DWOPER  Model 
CPU  Time 

5.10 

1.50  sec 

— 

5.11 

6 

400 

100 

6.57  sec 

256.64  sec 

5.12 

12 

100 

100 

0.49  sec 

— 

5.13 

12 

400 

100 

2.17  sec 

49.08  sec 

Double-Peaked  Hydrograph 


Figure 

°in 
(hours) 

5.14 

— 

5.15 

— 

5.16 

— 

5.17 

— 

L 
(miles) 

% 
(1000  cfs) 

100 

Nonlinear  State 

Space  Model 

CPU  Time 

DWOPER  Model 
CPU  Time 

100 

1.10  sec 

— 

400 

100 

4.36  sec 

77.25  sec 

400 

10 

4.34  sec 

77.25  sec 

400 

200 

4.35  sec 

77.25  sec 
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orders  of  magnitude  more  CPU  time  than  does  the  nonlinear  state  space 
model. 

Summary 

Results  from  the  nonlinear  state  space  and  DWOPER  models  have  been 
compared  for  large  inflow  transients  in  rectangular  and  triangular 
prismatic  channels.   Inflow  hydrographs  of  varying  durations,  and  a 
double-peaked  inflow  were  routed  through  channels  100  and  400  miles 
long.   The  results  presented  in  this  chapter  have  demonstrated  that: 

1.  Many  of  the  important  nonlinearities  in  the  full  Saint-Venant 
equations  are  emulated  by  the  nonlinear  state  space  model. 

2.  The  equations  for  computing  average  flow  in  a  reservoir 
generate  flow  levels  that  are  too  high  for  steeply  rising 
inflow  hydrographs.   This  leads  to  improper  state  space  model 
parameter  selection  and  degraded  model  performance. 

3.  The  nonlinear  state  space  model  approximates  well  the  outflow 
hydrograph  recession  of  the  DWOPER  model  for  all  inflow 
hydrographs  presented. 

4.  The  problem  in  Chapter  III  of  an  outflow  hydrograph  that  is 
too  symmetrical  is  solved  with  the  nonlinear  state  space 
model . 

5.  The  nonlinear  state  space  model  typically  preserves 
continuity  within  plus  or  minus  3  percent. 

6.  The  DWOPER  model  requires  approximately  2  orders  of  magnitude 
more  CPU  time  than  the  nonlinear  state  space  model. 

The  Saint-Venant  equations  have  been  our  exemplar  model  to 
illustrate  the  new  perspective  of  routing  in  hydrology.   A  general 
strategy  to  emulate  any  primary  model  is  the  subject  of  Chapter  VI. 


CHAPTER  VI 
HOW  TO  EMULATE  ANY  PRIMARY  MODEL 


Introduction 


Results  presented  in  Chapter  V  demonstrate  that  the  the  nonlinear 
state  space  routing  model  can  emulate  the  full  Saint-Venant  equations. 
This  is  only  one  application  of  the  general  approach  to  nonlinear 
modelling.   The  nonlinear  state  space  model's  structure  should  allow  it 
to  imitate  the  behavior  of  any  surface  runoff  or  channel  routing 
model.   This  chapter  describes,  in  detail,  the  steps  necessary  to 
emulate  any  primary  model.   Figure  6.1.  summarizes  those  steps. 

Select  a  Primary  Model 


The  obvious  first  step  is  to  choose  a  primary  model.   The  primary 
model  is  chosen  based  on  the  data  available  to  describe  the  physical 
system  and  the  computational  resources  committed  to  simulate  the 
catchment  or  channel. 

Calibrate  the  Primary  Model 


The  primary  model  must  be  calibrated  to  the  channel  or  catchment 
of  interest.   In  general,  the  parameters  of  the  nonlinear  state  space 
model  are  determined  based  on  simulation  with  the  primary  model.   The 
behavior  of  the  nonlinear  state  space  model  emulates  the  primary  model, 
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so  any  errors  in  calibration  of  the  primary  model  will  affect  the 
performance  of  the  nonlinear  state  space  model. 

Determine  the  Impulse  Response  Properties  of  the  Primary  Model 

Ideally  we  would  determine  the  weighting  function  that  the  primary 
model  uses  to  transform  input  signals  into  output  signals. 
Unfortunately,  it  is  physically  impossible  to  observe  the  weighting 
function  of  a  system.   As  described  in  Chapter  II,  for  a  linear 
stationary  system  the  weighting  function  is  the  time-reversed  image  of 
the  system  impulse  response.   The  premise  of  this  work  is  that  the 
mathematical  behavior  of  a  nonlinear  stationary  system  can  be 
represented  by  a  linear  model  whose  impulse  response  matches  that  of  the 
nonlinear  system  at  all  levels  of  flow.   The  impulse  response  of  a 
nonlinear  system  changes  with  flow  level;  if  it  did  not,  the  system 
would  be  linear.   To  emulate  the  changing  impulse  response  of  the 
nonlinear  system,  the  linear  model  must  be  nonstationary  (i.e.,  its 
parameters  and,  therefore  its  impulse  response,  must  vary  with  flow 
level) . 

Study  of  the  variations  in  the  mean  and  variance  of  the  impulse 
response  function  described  in  Chapter  IV  indicates  that  within  flow 
ranges  where  the  channel  shape  is  regular,  the  moments  of  the  impulse 
response  function  of  the  Saint-Venant  equations  change  as  simple 
algebraic  powers  of  the  kinematic  wave  parameter,  m.   Therefore,  the 
value  of  ra  can  be  found  from  the  variations  with  baseflow  level  of  the 
moments  of  the  system  impulse  response.   Equations  2.53  and  2.54  express 
the  mean  and  variance  of  the  system  impulse  response  function.   For  the 
mean  or  variance,  the  moment  of  the  impulse  response  function  is  equal 
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to  the  difference  of  the  moment  of  the  system  output  and  the  system 
input.   If  small  inflow  transients  (i.e.,  10  percent  of  baseflow)  are 
input  to  the  primary  model  at  several  baseflow  levels,  the  means  and 
variances  of  the  inflow  and  outflow  can  be  determined  for  each  baseflow 
level.   Transient  inflow  hydrographs  can  be  produced  with  the  Fortran 
program  listed  in  Appendix  B.   Inputs  to  this  program  are: 

1.  Q  ,  the  baseflow  level, 

2.  0p,  the  incremental  transient  flow  level, 

3.  Q  »  the  flow  level  about  which  the  equations  are  linearized, 

4.  a^n,  the  standard  deviation  of  the  hydrograph  produced, 

5.  Sn ,  B„  and  m,  which  define  the  channel  geometry,  and 

6.  the  time  step  and  total  period  for  which  the  hydrograph  will 
be  generated. 

In  general,  the  means  and  variances  of  the  impulse  response  will 
vary  with  flow  as  shown  in  Figures  4.3  through  4.6.   The  value  of  m  can 
be  found  from  the  slope  of  the  system  impulse  response  mean  (the  lag 
introduced  by  the  channel)  versus  discharge  relationship  in  log  space 
with  equation  4.7.   The  slope  of  the  lag  versus  flow  level  in  log  space 
may  change  with  flow  level  for  a  general  channel  shape.   Abrupt  changes 
in  channel  geometry  can  introduce  nonlinearities  not  accounted  for  with 
the  algebraic  power  relationship  of  impulse  response  properties  with 
flow.   In  fact,  for  a  trapezoidal  prismatic  channel,  the  value  of  m  must 
vary  with  flow.   For  shallow  flows,  a  trapezoidal  channel  behaves  like  a 
rectangular  channel  with  an  m  value  of  5/3.   As  the  flow  deepens,  the 
trapezoidal  channel  operates  more  like  a  triangular  channel  and  has  an  m 
value  of  4/3. 
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The  value  of  m  can  generally  be  expressed  In  a  piecewise  constant 
versus  flow  manner  as  in  Figure  6.2.   The  flow  ranges  over  which  m  is 
constant  can  be  determined  by  the  flow  levels  where  the  slope  is 
constant  in  the  log(lag)  versus  log(discharge)  relationship. 

For  the  particular  application  explored  in  Chapter  V,  the  value  of 
m  was  constant  for  all  flows  and  could  be  obtained  given  the  prismatic 
channel  shape.   For  any  prismatic  channel  shape  that  obeys  equation 
2.23,  the  value  of  m  can  be  found  from  Manning's  friction  formula.   In 
addition,  the  mean  and  variance  of  the  system  impulse  response  are 
expressed  as  functions  of  reference  flow  level  by  equations  2.56  and 
2.58. 

Estimate  a,  n  and  k  as  Functions  of  Flow  Level 

The  final  pieces  of  information  needed  for  the  nonlinear  state 
space  model  are  reference  values  of  the  model  parameters  ( a^  ,  r^  and  kg) 
at  any  flow  level  Q   in  each  of  the  constant  m  ranges.   The  variation  of 
parameters  with  flow  is  specified  by  the  value  of  m.   The  single 
reference  point  at  flow  level  Q   fixes  the  functions  for  a,  n  or  k.   The 
values  of  a  ,  a,  and  k.  can  be  found  from  the  impulse  response  mean  and 
variance  for  a  small  transient  input  on  a  baseflow,  Q  within  the  range 
of  flows  for  which  m  is  constant.   The  reference  parameter  values  are 
computed  with  equations  2.78  through  2.80. 

Parameters  of  the  Nonlinear  State  Space  Model 

The  parameters  of  the  nonlinear  state  space  model  are  now 
specified.   Given  m  and  a  f    n     and  kQ  at  flow  level  CL  for  each  constant 
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m  flow  range,  any  inflow  transient  can  be  modelled  with  the  Fortran 
program  listed  in  Appendix  C.   The  behavior  of  the  nonlinear  state  space 
model  mathematically  emulates  the  primary  model. 


CHAPTER  VII 
RELATED  PREVIOUS  WORK 


Introduction 


Many  other  authors  have  approached  the  problem  of  routing  surface 
runoff  or  channel  flow.   This  chapter  describes  pertinent  previous  work 
in  the  field  of  hydrology.   Surface  runoff  and  channel  routing  models 
related  to  the  nonlinear  state  space  model  by  either  their  common 
cascade  of  reservoirs  structure  or  approach  of  modelling  nonlinearities 
by  varying  parameter  values  are  described  in  the  next  section.   Observed 
flow  level  versus  lag  time  data  which  corroborates  the  theoretical 
relationships  developed  in  Chapter  IV  are  presented.   In  addition,  the 
nonlinear  state  space  and  Muskingum-Cunge  models  are  compared  with  the 
DWOPER  model  for  a  400  mile  long  rectangular  channel. 

Related  Routing  Models 

Early  work  with  models  constructed  as  cascades  of  reservoirs  is 
attributed  to  Nash  (L957)  and  to  Kalinin  and  Miljukov  (1958).   Many 
authors  followed  this  work  in  attempts  to  model  a  catchment 's 
instantaneous  unit  hydrograph  (IUH)  (Rao,  Oelleur  and  Sarma  1972;  Sauer 
1973;  Pedersen,  Peters  and  Helweg  1980).   The  IUH  is  the  common  term  in 
the  field  of  hydrology  for  a  system's  impulse  response  function.   More 
recent  work  with  cascades  of  reservoirs  has  dealt  with  the  development 
of  a  state  space  model  of  the  Nash  cascade  (Szollosi-Nagy  1981).   The 
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above-mentioned  works  all  deal  with  linear  reservoir  models  that  are 
strictly  applicable  only  over  a  small  flow  range  around  the  reference 
discharge. 

Amorocho  and  Brandstetter  (1971),  and  Napiorkowski  and 
Strupczewski  (1981)  worked  with  integral  forms  of  the  open  channel  flow 
equations  to  parameterize  the  properties  observed  in  the  IUH  when 
routing  unsteady  flow.   The  resultant  mathematical  forms  were  cumbersome 
and  not  amenable  to  state  space  formulation.   Keefer  and  McQuivey  (1974) 
took  the  approach  of  splitting  the  flow  range  into  sections  and  forming 
a  linear  model  in  each  range.   Calibration  of  the  multiple  linearization 
model  was  complicated  because  an  impulse  response  function  had  to  be 
found  for  each  flow  range.   They  recognized  that  the  IUH  varies  with 
flow,  but  did  not  study  the  form  of  that  variation. 

Other  authors  developed  simple  nonlinear  models  of  surface  runoff 
or  channel  flow  (Singh  1964;  Dooge  1967;  Prasad  1967;  Reed,  Johnson  and 
Firth  1975;  Singh  and  Buapeng  1981).   These  studies  of  the  IUH  generally 
expressed  the  impulse  response  properties  in  terms  of  basin 
characteristics  through  complex  regression  formulas.   Reed,  Johnson  and 
Firth  comment  that  further  study  to  determine  the  form  of  the 
relationship  between  physical  factors  and  variable  lag  is  needed. 

Several  authors  let  model  parameters  vary  with  time  rather  than 
flow  level.   Mandeville  and  O'Donnell  (1973)  developed  time  varying 
expressions  for  linear  reservoir  and  linear  channel  equations.   The  time 
varying  functions  for  reservoir  and  channel  parameters  are  fit  to 
observed  data.   A  model  in  which  the  parameters  change  with  the  season 
of  the  year  was  developed  by  Chiu  and  Bittler  (1967).   They  attempted  to 
model  the  observed  nonlinear  phenomena  by  looking  for  cycles  in  flow 
patterns  from  year  to  year. 
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Mein,  Laurenson  and  McMahon  (1974)  presented  a  nonlinear  storage 
reservoir  model  with  the  proper  structure  to  simulate  the  nonlinearities 
in  surface  runoff  or  channel  flow.   The  variation  with  flow  of  the 
nonlinear  properties  is  empirically  fit  to  observed  data.   Their  model 
cannot  easily  be  put  in  state  space  form. 

Schaake  (1965,  pp.  76-92)  recognized  that  the  IUH  derived  from  the 
rising  portion  of  a  large  surface  runoff  event  was  different  than  the 
IUH  from  the  recession,  and  that  the  impulse  response  function  varied 
with  the  magnitude  of  an  event.   He  concluded  that  the  differences  were 
caused  by  nonlinear  properties  of  flow  routing  for  large  transients. 
Because  of  Schaake 's  analysis,  Valdes,  Fiallo  and  Rodriguez-Iturbe 
(1979)  studied  incremental  IUHs  and  concluded  that  there  were  not 
significant  differences  in  the  IUHs  computed  for  the  rising  and  falling 
limbs  of  hydrographs  for  small  inflow  transients. 

None  of  the  models  discussed  above  has  the  important  combination 
of  properties  found  in  the  nonlinear  state  space  model  developed  in  this 
dissertation.   The  analysis  conducted  in  the  current  work  has  recognized 

that 

1.  the  variations  with  flow  of  model  parameters  can  be  described 
functionally  by  relating  the  parameters  to  the  mean  and 
variance  of  the  incremental  impulse  response  for  various  flow 
levels,  and 

2.  a  state  space  structure  is  invaluable  for  second-order 
analysis  by  means  of  estimation  theory. 
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Observed  Variations  with  Flow  Level  of  the  System  Lag  Time 

The  variation  of  lag  with  flow  level  has  been  known  in  hydrology 
for  some  time.   Laurenson  (1964)  empirically  derived  a  relationship 
between  lag  and  flow  level  by  analyzing  hydrographs  from  the  South  Creek. 
Experimental  Catchment  near  Sydney,  New  South  Wales,  for  1956  through 
1960.   Laurenson  plotted  the  lag  versus  mean  flow  values  on  linear,  then 
semi-log,  and  finally  log-log  axes.   He  found  that  the  data  essentially 
describe  a  straight  line  in  log-space.   The  data  he  used  and  the 
regression  equation  he  fit  to  the  data  are  shown  in  Figure  7.1.   The 
slope  of  the  regression  line  he  fit  in  log-space  is  -3.70.   From 
equation  4.22  we  obtain  a  value  of  m  =  1.37  for  a  slope  of  -3.70.   This 
value  of  m  falls  between  those  of  a  rectangular  (m  =  1.67)  and  a 
triangular  (m  =  1.25)  channel. 

Askew  (1970)  continued  on  the  same  path  as  Laurenson,  recognizing 
the  relationship  between  lag  time  and  flow  level.   He  fit  slightly 
different  data  than  did  Laurenson  for  South  Creek,  and  determined 
regression  equations  for  lag  versus  mean  flow  for  four  other  catchments 
in  Australia.   The  regression  equation  exponents  and  the  associated  m 
values  found  with  equation  4.22  are  presented  in  Table  7.1.   All  the 
values  for  m  fall  close  to  or  within  the  1.25  to  1.67  range  expected 
theoretically  for  triangular  and  rectangular  cross  section  shapes. 

Comparison  with  the  Muskingum-Cunge  Model 

The  Muskingum-Cunge  routing  model  is  presented  in  Appendix  A  as 
equations  A. 7  through  A. 9.   This  model  is  developed  from  a  storage 
routing  equation  using  kinematic  wave  theory  (Cunge  1969).   The 
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Figure  7.1  Lag  Versus  Mean  Discharge  Relationship, 
South  Creek  Experimental  Catchment, 
University  of  New  South  Wales 
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Table  7.1   Values  of  m  for  Several  Catchments 


Slope 

from  Askew 

Catchment 

Regress 

ion  Equa 
-4.149 

tion 

m 

South  Creek 

1.32 

Hacking   River 

-5.181 

1.24 

Eastern  Creek 

-3.802 

1.36 

Cawleys   Creek 

-5.263 

1.23 

Research  Creek 

-3.049 

1.49 
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Muskingum-Cunge  model  Is  the  most  versatile  and  physically  relavent  of 
the  storage  routing  models  (Fread  1982). 

The  nonlinear  state  space  and  Muskingum-Cunge  models  are  compared 
with  the  DWOPER  model  in  Figure  7.2.   A  single-peaked  inflow  transient 
was  routed  through  a  400  mile  long  rectangular  channel  with  each  of  the 
three  models.   Parameters  of  the  nonlinear  state  space  model  are  the 
same  as  in  Chapter  V,  for  400  mile  long  rectangular  channels. 
Parameters  of  the  Muskingum-Cunge  model  are  derived  from  channel 
properties  as  shown  in  equations  A. 8  and  A. 9. 

The  Muskingum-Cunge  model  matches,  in  general,  the  shape  of  the 
outflow  hydrograph  produced  by  the  DWOPER  model.   However,  results  with 
the  Muskingum-Cunge  model  peak  slightly  late  and  high,  and  remain  above 
the  DWOPER  model  outflow  throughout  the  recession.   The  nonlinear  state 
space  model  is  early  on  the  rising  limb,  but  the  timing  of  the  peak  flow 
and  the  recession  match  the  DWOPER  model  excellently. 

The  ratio  of  CPU  time  required  by  the  Muskingum-Cunge  model  to 
that  of  the  nonlinear  state  space  model  is  approximately  25:1.   The 
distance  step  used  in  the  numerical  solution  of  the  Muskingum-Cunge 
equations  was  5  miles;  implying  that  the  Muskingum-Cunge  model  used  27 
subreaches  to  simulate  the  example  channel.   The  average  number  of 
reservoirs  needed  by  the  nonlinear  state  space  model  was  8.057.   This 
ratio  of  the  number  of  distance  steps  to  reservoirs  of  3:1  is  consistant 
with  results  presented  in  Chapter  V,  and  is  undoubtedly  a  contributing 
factor  to  the  additional  CPU  time  needed  by  the  Muskingum-Cunge  model. 
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CHAPTER  VIII 
CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FUTURE  WORK 


Conclusions 

The  work  presented  herein  consists  of  a  new  perspective  of  routing 
in  hydrology  and  its  application.   The  new  perspective  was  developed  in 
Chapter  I.   A  particular  primary  model  was  then  chosen  to  demonstrate 
the  applicability  of  the  new  perspective.   The  exemplar  primary  model 
was  the  set  of  equations  of  unsteady  flow  in  open  channels.   General 
conclusions  about  the  new  perspective  are: 

1.  The  essence  of  this  new  perspective  of  routing  is  its  two-step 
approach  to  the  problems  encountered  with  existing  routing 
technology.   First,  the  physics  of  a  system  are  modelled  with 
a  primary  model;  then  the  mathematics  of  the  primary  model  are 
emulated  by  a  state  space  model.   The  state  space  model  is  a 
model  of  the  primary  model. 

2.  Stationary  nonlinear  behavior  of  a  primary  model  can  be 
emulated  by  a  nonstationary  linear  model  having  a  state  space 
formulation. 

3.  The  nonlinearities  present  in  the  primary  model  can  be 
captured  by  matching,  at  various  reference  flow  levels,  the 
incremental  impulse  response  properties  of  the  primary  model 
with  a  linear  state  space  model  structure. 

4.  Using  a  primary  model  to  account  for  the  observed  physics,  and 
mathematically  emulating  the  primary  model's  behavior  with  a 
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state  space  model  whose  parameters  vary  with  flow  level  is 
appropriate  for  modelling  surface  runoff  and  channel  routing 
phenomena  for  downstream  flow  in  prismatic  channels. 

In  addition,  several  conclusions  derived  from  the  application  of  the  new 

perspective  to  the  Saint-Venant  equations  are: 

1.  The  impulse  response  function  of  the  linearized  Saint-Venant 
equations  in  a  prismatic  channel  has  the  form  of  an  inverse 
Gaussian  pdf. 

2.  For  small  (10  percent  of  baseflow)  inflow  transients,  the 
inverse  Gaussian  pdf  models  well  the  behavior  of  the  full 
equations  of  unsteady  flow. 

3.  The  parameters  of  the  inverse  Gaussian  pdf  can  be  computed 
from  the  moments  of  the  channel  impulse  response  function, 
which  is  determined  by  simulating  with  the  primary  model  for 
small  inflow  transients. 

4.  For  a  prismatic  channel,  the  inverse  Gaussian  pdf  parameters 
can  be  found  directly  as  functions  of  the  channel  geometry. 

5.  Parameters  of  the  inverse  Gaussian  pdf  are  related  in  a  very 
simple  way  to  parameters  of  the  gamma  pdf. 

6.  A  state  space  model,  based  on  the  concept  of  a  cascade  of 
reservoirs,  has  parameters  identical  to  the  3-parameter  gamma 

pdf. 

7.  A  linear  state  space  model  (i.e.,  one  whose  parameters  are 
constant  with  flow  level)  is  an  adequate  surrogate  for  the 
full  Saint-Venant  equations  for  small  inflow  transients  in  a 
prismatic  channel,  but  its  behavior  deterioriates  when  large 
transients  are  routed. 
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8.  The  mean  and  variance  of  the  impulse  response  function  of  the 
equations  of  unsteady  flow  vary  as  simple  algebraic  powers  of 
flow  level;  the  expressions  for  the  exponents  are  functions  of 
the  channel  geometry. 

9.  A  nonlinear  state  space  model,  whose  parameters  vary  to  match 
the  changes  with  flow  level  of  the  impulse  response  function's 
mean  and  variance,  is  an  excellent  surrogate  for  the  NWS 
DWOPER  model,  which  is  a  numerical  solution  of  the  full 
Saint-Venant  equations. 

10.  The  DWOPER  model,  which  numerically  solves  the  Saint-Venant 
equations  in  a  very  efficient  manner,  requires  approximately  2 
orders  of  magnitude  more  CPU  time  on  a  PRIME  750  computer  than 
does  the  nonlinear  state  space  model. 

11.  The  ability  of  the  state  space  model  to  emulate  the  full 
Saint-Venant  equations  deteriorates  for  very  abrupt  inflow 
hydrographs,  such  as  from  a  sudden  dam  failure,  or  a  small, 
steep-sloped  catchment. 

12.  The  current  form  of  the  nonlinear  state  space  model  is  not 
suitable  for  modelling  channels  with  lateral  inflow  or 
tributaries;  however,  this  limitation  should  be  easily 
overcome . 

13.  Backwater  conditions  and  upstream  wave  movement  are  not 
accounted  for  because  of  the  underlying  structure  of  the  state 
space  model,  which  is  based  on  the  diffusion  equation. 
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Recommendations  for  Future  Work 

Some  of  the  inadequacies  of  current  routing  technology  listed  in 
Chapter  I  are: 

1.  it  produces  no  objective  estimate  of  uncertainty, 

2.  repetitive  computations  such  as  in  flood  frequency  analyses 
can  be  costly  for  fully  nonlinear  models,  and 

3.  error  properties  of  models  and  data  are  not  expressly 
considered. 

The  current  work  directly  addresses  only  item  2  above.   The  fully 
nonlinear  state  space  model  developed  in  this  work  is  several  orders  of 
magnitude  more  CPU  time  efficient  than  a  numerical  solution  of  the 
Saint-Venant  equations.   Although  the  state  space  model  has  a  structure 
which  is  amenable  to  the  application  of  estimation  theory,  and  could 
produce  an  objective  estimate  of  uncertainty  and  consider  error 
properties,  that  step  remains  for  future  work.   Several  intermediate 
tasks  are  required  to  formulate  a  filter  theory  model  from  the  nonlinear 
state  space  model.   These  tasks  include: 

1.  developing  a  stage  versus  discharge  (rating  curve)  model  which 
includes  a  stochastic  component  to  account  for  shifts  in  the 
rating  curve, 

2.  determining  a  technique  to  specify  values  for  the  filter 
parameters,  and 

3.  addressing  the  theoretical  problem  of  filtering  when  the 
number  of  model  states  varies  with  time. 

The  nonlinear  state  space  model  routes  discharge  on  a  catchment  or 
in  a  channel.   Observations  of  hydrographs  are  most  often  made  in  terms 
of  stage  or  depth  of  flow.   The  relationship  between  stage  and  discharge 
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is  not  a  simple  one  (Fread  1973).   For  a  given  stage,  the  discharge  is 
higher  as  the  rising  limb  of  a  hydrograph  passes  a  given  point  than 
during  the  recession.   This  is  because,  during  the  rising  limb,  the 
water  level  upstream  is  relatively  higher  than  downstream,  so  the  forces 
moving  the  water  along  the  channel  are  greater.   During  a  recession, 
with  the  wave  crest  now  downstream,  the  water  surface  slope  is  less  and, 
therefore,  a  smaller  force  is  acting  on  the  water.   In  addition  to  the 
variations  in  rating  curves  due  to  water  surface  slope,  channel  bed  form 
changes  can  have  dramatic  effects,  especially  in  sand  bed  rivers 
(Simons,  Stevens  and  Duke  1973).   A  suitable  rating  curve  model 
structure  compatible  with  the  nonlinear  state  space  routing  model  and 
yet  accounting  for  the  uncertainties  inherent  in  the  stage-discharge 
relationship  must  be  developed. 

Kalman  filtering  is  a  technique  that  combines  the  information  in  a 
model  and  measurements  of  a  process  to  provide  a  better  estimate  of  the 
current  conditions  of  the  process  than  either  the  model  or  measurements 
could  give  alone  (Kalman  1960).   A  Kalman  filter  algorithm  can  be 
implemented  as  equation  1.1 


and 


X    =  F    'X   +  G    «U   , 
-t+1    -X,t  -t   -X,t  -t+1 


Y   =  H  *X  +  V  (8.1) 

_t   _t  _t   -t 


where  X,  F,  G  and  U_ are  defined  as  in  equation  1.1, 


e{uc}   =  0 


E{V        =   0 
L— tJ 


EivWi  =v6i 
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in  which  Y^  is  the  vector  of  observations  at  time  t, 

Hf  is  the  measurement  matrix, 

Vj.  is  the  vector  of  measurement  errors, 

q,.  is  system  noise  covariance  matrix, 

Rj.  is  the  measurement  noise  covariance  matrix, 

T  indicates  a  matrix  transpose, 

E  is  the  expectation  operator,  and 

6 .  is  the  Dirac  delta  function  which  is  equal  to  1  when 

i  =  0  and  otherwise  is  zero  (Jazwinski  1970,  p.  269). 

Equation  1.1  is  the  structure  with  which  the  nonlinear  state  space  model 
was  developed.   In  a  Kalman  filtering  algorithm  it  is  called  the  system 
equation  and  is  responsible  for  propagating  the  system  dynamics  in 
time.   Equation  8.1  is  called  the  measurement  equation  and  is  used  to 
incorporate  information  from  measurements  of  the  state  vector,  X_,  into 
the  filter  algorithm.   The  filter  parameters  that  must  be  specified  are 
the  matrices  0_  and  R_.   Determining  values  for  Q_  and  R_will  require  a 
major  portion  of  the  effort  to  form  an  estimation  theory  model  of 
surface  runoff  and  channel  routing  with  the  nonlinear  state  space  model. 

A  theoretical  issue  that  must  be  addressed  before  the  nonlinear 
state  space  model  can  be  used  in  a  Kalman  filter  algorithm  is  how  to 
account  for  a  varying  number  of  states.   The  technique  developed  in 
Chapter  IV  to  interpolate  the  outflows  from  each  reservoir  of  letting  n 
vary  works  well  when  one  is  simulating  with  equation  1.1  alone.   When, 
however,  the  Kalman  filter  is  processing,  it  propagates  not  only  the 
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state  vector,  but  also  the  covariance  matrix  of  the  error  in  X.     The 
system  error  covariance  matrix  has  a  dimension  of  n  by  n,  with  elements 
that  cannot  simply  be  averaged,  as  is  done  with  the  state  vector.   The 
first  step  in  attacking  this  problem  should  be  to  explore  the  literature 
of  other  physical  science  disciplines  where  estimation  theory  has  been 
extensively  employed. 

The  nonlinear  state  space  model  developed  in  this  dissertation  is 
based  on  a  continuous  function  analysis.   The  approach  taken  was  to  work 
with  the  gamma  pdf,  a  continuous  function.   Differential  equations  that 
match  the  gamma  pdf  were  developed  as  the  nonlinear  state  space  model 
structure.   This  restricted  us  to  small  time  steps  in  the  actual  model 
implementation  since  we  were  trying  to  emulate  a  continuous  function. 
An  alternative  approach  is  to  work  with  a  discrete  function  for  the 
incremental  impulse  response  function  (i.e.,  the  pulse  response 
function).   The  general  strategy  would  be  to  find  a  discrete  function 
with  the  analytical  form  of  the  pulse  response,  and  then  equate  the 
function  and  pulse  response  by  matching  their  moments.   Time  steps  for 
the  model  derived  from  a  discrete  analysis  would  be  limited  only  by  the 
frequency  properties  of  the  inflow  hydrograph,  not  by  numerical  aspects 
of  approximating  a  continuous  function. 

By  specifying  the  changes  with  flow  level  of  the  impulse  response 
function  in  terms  of  its  mean  and  variance,  we  have  limited  the  shape  of 
impulse  response  functions  that  can  be  approximated  by  the  nonlinear 
state  space  model  to  unimodal  ones.   A  different  parameterization  of  the 
impulse  response  function  could  allow  multi-peaked  system  responses  to 
be  simulated  with  a  nonstationary  linear  model  structure.   An  analysis 
in  the  spectral  domain  may  be  the  appropriate  approach  for  this  problem. 
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Surface  runoff  and  channel  routing  comprise  only  a  portion  of  the 
hydrologic  cycle  (Eagleson  1970,  pp.  5-8).   Subsurface  flow  routing  is 
another  critical  path  in  the  constant  motion  of  water.   The  highly 
nonlinear  nature  of  subsurface  flow  has  been  modelled  with  2-  and 
3-dimensional  versions  of  the  equations  of  motion  (Eagleson  1970,  pp. 
271-282).   Groundwater  flow  is  modelled  with  a  conceptual  structure  in 
the  NWS  Sacramento  Soil  Moisture  Accounting  (SMA)  model  (Burnash,  Ferral 
and  McGuire  1973,  pp.  11-12).  Models  in  state  space  form  have  been 
developed  for  the  NWS  SMA  model  (Kitanidis  and  Bras  1978,  Goldstein  and 
Larimore  1980).   These  models  retain  the  inherent  structure  of  the 
original  NWS  SMA  model  in  a  set  of  simultaneous  differential  equations 
used  as  the  system  equation,  1.1.   It  may  be  possible  to  extend  the  new 
perspective  of  routing  to  account  for  the  highly  nonlinear  structure  of 
groundwater  flow.   The  approach  is  not  clear,  however,  since  the 
nonlinear  properties  do  not  seem  to  be  simple  functions  of  flow  level, 
but  instead  are  related  to  the  interactions  of  many  components. 

A  final  area  for  future  work  is,  of  course,  to  use  the  nonlinear 
state  space  routing  model  on  real  catchments  and  channels.  Work,  with 
real  data  was  eschewed  in  this  dissertation  because  the  intent  of  the 
new  perspective  is  to  use  the  state  space  model  as  a  surrogate  for  a 
model  of  the  catchment  or  channel  physics.   The  objective  of  using  an 
exemplar  primary  model  was  to  develop  the  theoretical  basis  for  relating 
the  primary  and  state  space  models.   The  strategy  to  employ  the  new 
perspective  has  been  presented  in  Chapter  VI.   This  strategy  and  the 
Fortran  programs  listed  in  Appendices  B  and  C  must  be  tested  for  actual 
catchments  and  channel  reaches. 


APPENDIX  A 
A  REVIEW  OF  SOME  HYDROLOGIC  ROUTING  MODELS 


The  complete  Saint-Venant  equations  of  gradually  varied,  unsteady 
flow  in  an  open  channel  are 

B.8y  +  MAv)_  (A#1) 

3 1      dx 

and 

1   9v        3y        v  3v  q      r  \  tk    0\ 

gTFTxgTx  f  0         g«A  v    x  ; 

where        t    =   time, 

x  =  distance  along  the  channel, 

y  =  depth  of  flow, 

v  =  average  velocity  of  the  flow, 

A  =  cross  sectional  area  of  the  flow, 

B  =  surface  width  of  the  flow, 

g  =  the  acceleration  of  gravity, 

q  =  lateral  inflow  per  unit  length, 

u   =  the  x  component  of  the  velocity  of  the  lateral  inflow, 

Sf  =  the  friction  slope,  and 

S„  =  the  channel  bottom  slope  (Henderson  1966,  p.  287). 
The  conservation  of  mass  is  described  by  the  continuity  equation,  A.l. 
The  conservation  of  momentum  is  expressed  by  the  equation  of  motion, 
A. 2.   The  terms  of  equation  A. 2,  from  left  to  right,  are  measures  of  the 
local  acceleration,  the  pressure  force,  the  convective  acceleration,  the 
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friction  force,  gravity  and  the  acceleration  of  the  lateral  inflow.   The 
quasi-linear,  hyperbolic  system  of  equations  defined  by  equations  A.l 
and  A. 2  has  no  known  analytical  solution.   Numerical  solutions  are 
possible  when  initial  and  boundary  conditions  are  specified. 

All  channel  routing  models  obey  the  continuity  equation,  A.l,  but 
may  be  categorized  based  on  which  terms  of  equation  A. 2  they  consider. 
Certain  routing  models  do  not  consider  any  terms  of  equation  A. 2.   They 
are  generally  called  storage  routing  models  and  are  based  on  the 
continuity  equation  alone.   Equation  A.l  is  usually  rewritten  as 

I-O-J*  (A.3) 

At 

where   I  is  the  channel  inflow, 

0  is  the  channel  outflow,  and 

AS  is  the  change  in  storage  in  the  channel  during  a  At  time 

interval. 
The  storage,  inflow  and  outflow  may  be  related  by 

S  =  K'[x»I  +  (1-X)*0]  (A. 4) 

where   K  and  X  are  model  parameters. 

Reservoir  routing  models,  such  as  those  developed  by  Puis  (1928) 

and  Goodrich  (1931)  set  X  =  0  in  equation  A. 4.   Storage  depends  only  on 

outflow,  and  equation  A.3  can  be  expressed  in  centered  finite  difference 

form  as 

I+I  0+0  S-S 

_J L  -  _J 2.  =  _2 L  (A. 5) 

2  2  At 

where  the  subscripts  1  and  2  indicate  the  previous  and  current  time 
steps,  respectively.   Equation  A. 5  can  be  solved  for  the  unknowns  at  the 
current  time  step  as 
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2«S  2»S 

-r-?-  +0  =1  +1  +  -r-L  ~   0  (A. 6) 

At     2     1     2    At     1 

The  outflow  0,  can  be  found  given  an  S-0  relationship  from  observed 
data. 

The  well-known  Muskingum  model  can  be  found  by  substituting 
equation  A. 4  into  equation  A. 5  to  yield 

0   =  C  •  I  +  C  «I   +C«0  (A. 7) 

2     0   2     11     2   1 

where   CQ  =  -(K>  X  -  At/2)/C  , 

C:  =  (K«X  +  At/2)/C3  , 

C2  =  (K  -  K'X  -  At/2)/C3,  and 

C3  =  K  -  K'X  +  At/2  (Chow  1964,  p.  25-40). 

In  the  original  Muskingum  model  as  developed  by  McCarthy  (1938),  the 
parmaters  K  and  X  were  determined  from  observed  inflow-outflow  data. 
Cunge  (1969)  developed  the  Muskingum  model  from  kinematic  wave 
theory  with  a  single-valued  stage-discharge  relationship.   In  the 
Muskingum-Cunge  model  the  parameters  K  and  X  of  equation  A. 7  are 
specified  as 

K  =  Ax/c  (A. 8) 

and 

X  =  L[\   -  0  /(B  «c'S  -Ax)]  (A. 9) 

2       0    0     0 

where   c  =  — ,  is  the  kinematic  wave  speed  at  reference  discharge  0  , 
dA  v 

Ax  =  the  reach  length, 

S      =   the  channel  bottom  slope,  and 

B   =  the  channel  width  at  QQ . 
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Kalinin  and  Miljukov  (1958)  developed  another  model  related  to  the 
Muskingum  model,  which  is  expressed  as 

°2  =  °1  +  ^1  "  °1^K1  +  ^2  "  Xl^K2  (A*10) 

where   K  =  1  -  e-<c*At>/Ax, 

1 

K  =  1  -  K  •Ax/Cc'At),  and 
2         1 

Ax  =  Q  /S  •  (Ah/AQ) 
0   0 

in  which  Ah/AQ  is  the  slope  of  the  stage-discharge  relationship. 

The  Muskingum  model  is  equivalent  to  equation  A. 10  If  K  =  Ax/c  and  X  =  0 

(Fread  1982). 

In  the  Lag  and  K  model  (Linsely,  Kohler  and  Paulhus  1949,  pp.  230- 

232)  the  inflow  is  first  lagged,  or  delayed,  and  then  attenuated  with 

the  expression 

0   =  [i   +1   -  0  «(1  -  2-K/At)]/(l  +  2'K/At)  (A. 11) 

2      1     2     1 

This  model  is  linear  if  the  lag  and  K  factors  are  constant.   The  lag  and 
K  may  vary  with  flow  level  to  give  an  empirical  nonlinear  model. 

Complete  dynamic  models  retain  all  the  terms  of  equation  A. 2. 
They  provide  the  most  accurate  models  available,  hut  must  still  be 
considered  somewhat  empirical  in  nature,  since  the  effects  of  varying 
channel  geometry  are  lumped  into  an  empirical  roughness  coefficient 
(Weinmann  and  Laurenson  1979).   In  most  practical  applications  the  value 
of  Sn  is  several  orders  of  magniture  larger  than  the  acceleration  terms 
of  equation  A.l  (Henderson  1966,  p.  364).   A  reasonable  approach  to 
simplify  the  full  Saint-Venant  equations  is  to  neglect  the  acceleration 
terms  of  equation  A. 2.   Various  models  which  are  simplified  in  this  way 
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are  called  approximate  dynamic,  diffusion  analogy  or  kinematic  wave 
models  (Weinmann  and  Laurenson  1979).   Some  of  the  similarities  among 
the  approximate  routing  models  can  be  seen  by  rearranging  the  momentum 
equation,  A. 2  into  a  looped  rating  curve  formula.   This  is  accomplished 
by  beginning  with  a  flow  resistance  formula  in  general  form 


Q  =  C-A-R1"-/"^ 


(A.12) 


where   C  =  an  empirical  resistance  coefficient, 

R  =  the  hydraulic  radius,  and 

m  =  and  empirical  coefficient. 
For  steady,  uniform  flow,  Qu,  we  can  substitute  SQ  for  Sf  to  give 


Q  =  C«A»Rm»/S 


(A. 13) 


By  substituting  into  equation  A. 13  for  O A* Rm  from  equation  A. 12,  for  Sf 
from  equation  A. 2,  and  ignoring  the  lateral  inflow  term  we  can  derive  a 
general  looped  rating  curve  formula 


Q  =  Q  {  1  - 


1  3y    v  3v  _   1   3v;  1/2 
S^ST  "  S  g  57   S-^ 


Ht 


(A. 14) 


1   ° 

kinematic  wave        ! 
(i.e.,  steady,  uniform  flow) 

diffusion  analogy 

0 

0 

steady,  nonuniform  flow 

complete  dynamic  wave  model 

(i.e.,  unsteady,  nonuniform  flow) 

The  kinematic  and  diffusion  analogy  approximations  to  the  full 
Saint-Venant  equations  are  found  by  including  the  indicated  terms  of 
equation  A. 14  (Henderson  1966,  p.  287,  Weinmann  and  Laurenson  1979). 

Approximate  dynamic  models  are  based  on  the  continuity  equation 
expressed  as 
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»♦»§.,  (A.i5) 

and  the  specified  part  of  equation  A. 14.   Combining  equation  A. 15  and 
the  diffusion  analogy  terms  of  equation  A. 14  produces  the  equation  for  a 
convective  diffusion  model 

80     3Q     82Q  /»  wn 

*rt  +  C'r-  =  D« — -  +   c«q  (A. 16) 

at   H-   3x2 

where  for  regular  channels 

c  =  tt'-t^  ,  and 
B  dy  ' 


D  = 


2*B»S 

0 


The  parameters  c  and  D  control  the  convective  and  attenuation 
properties,  respectively.   D  is  called  the  diffusion  coefficient. 
Equation  A. 16,  without  the  lateral  inflow  term,  approximates  equation 
2.35  for  small  Froude  numbers.   Equation  A. 16  was  originally  proposed  by 
Hayami  (1951)  with  constant  parameters  c  and  D.   This  is  essentially  the 
equation  solved  by  the  linear  state  space  model  developed  in  this 
dissertation.   By  allowing  the  parameters  c  and  D  to  vary,  equation  A. 16 
describes  a  nonlinear  model  (Price  1973).   This  is  the  concept  behind 
the  nonlinear  state  space  model.   Both  Cunge  (1969)  and  Koussis  (1976) 
have  demonstrated  that  with  proper  choice  of  parameters,  the  Muskingum 
model  is  a  second  order  approximation  of  the  diffusion  equation. 

For  a  kinematic  wave  model  the  discharge  is  always  a  single-valued 
function  of  depth  (i.e.,  the  steady,  uniform  discharge)  because  from 
equation  A. 14 

Q  =  Q  (A. 17) 
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The  kinematic  wave  model  does  not  account  for  wave  attenuation.   The 
kinematic  wave  equation  can  he  derived  by  letting  D  =  0  in  equation  A. 16 
and  rearranging  to  give 

where   c  is  now  called  the  kinematic  wave  speed. 

If  c  is  constant,  equation  A. 18  describes  a  linear  model.   A  nonlinear 

kinematic  wave  model  arises  if  c  varies. 

Linearization  of  equation  A. 2  about  a  reference  flow  level  leads 
to  several  other  routing  models.   A  form  of  the  well-known  telegraph 
equation  can  be  obtained  by  linearizing  the  friction  term,  neglecting 
the  convective  acceleration,  and  assuming  zero  bottom  slope  and  a 
rectangular  channel,  as 

92v   32v     _  3v  /.  ,Q>. 

g.y. =  __  +  g.C  -J-  (A. 19) 

3X2    3t2      0  W 
where   CQ  ±s   a  constant  related  to  the  linearized  friction  term 
(Fread  1982). 
Lighthill  and  Whitham  (1955)  and  Harley  (1966)  developed  linear 
models  of  the  full  Saint-Venant  equations.   The  complete  linear  equation 
obtained  by  Harley  is 


32q    ^     92q     d2 


y3  _  v2).^JL-  2-v  -3-— -  

0     QJ    9x2       0  ^'3t    3t2 

(A. 20) 


u  0 


where   q  is  a  unit  discharge  in  a  unit  width  channel, 
S„  is  the  bottom  slope,  and 
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the  Saint-Venant  equations  have  been  linearized  about  a 

reference  velocity,  v  =  qQ /y0  (Harley  1966). 
Equation  A. 20  is  a  particular  case  of  equation  2.33  for  a  rectangular 
channel  with  m  =  3/2.   Because  Harley  used  a  Chezy  friction  factor  in 
the  derivation  of  equation  A. 20,  the  value  of  m  =  3/2  is  appropriate  for 
a  rectangular  channel.   The  impulse  response  function  of  equation  A. 20 
is 


h(x,t)  =  e  P'X'6(t-x/C  )  + 


n-fx/C     -  x/C  >eS'x  r#t«l{2'TTm}/m 

^12 


(A. 21) 


where   C  =  v  +  /g*y   , 

1  0        0 

C  =  v  -  /g«y   , 

2  0        0 


F  =  v  //g«y  , 

0      0 

P  =  S  «(2-F)/[2-y  -(f2+f)]  , 
0  0 

r  =  S  -v  -(2+F2)/(2-y  •  F2  J  , 
0   0  0 

s  =  SQ/(2.yo]  , 

m  =/(t-x/C  )-(t-x/C  )  ' 
1         2 

l{'}  is  a  first  order  Bessel  function  of  the  first  kind,  and 

fi(«)  is  the  Dirac  delta  function. 
Equation  A. 20  is  related  to  the  diffusion  analogy  model  of  Hayami  (Chow 
1959,  pp.  601-604).   As  with  any  linear  model,  the  response  of  the 
complete  linearized  model  is  very  dependent  on  the  reference  flow  level 
(Bravo,  Harley,  Perkins  and  Eagleson  1970,  pp.  35-39). 


APPENDIX  B 

A  FORTRAN  PROGRAM  TO  GENERATE  HYDROGRAPHS  WITH 

INVERSE  GAUSS UN  PDF  SHAPES  IN  A  PRISMATIC  CHANNEL 


c 

c 

C  THIS  PROGRAM  GENERATES  AN  INVERSE  GAUSSIAN  SHAPE  HYDROGRAPH 

C  FOR  RECTANGULAR  (EM-5/3)  OR  TRIANGULAR  (EM*4/3)  CHANNELS  FOR 

C  SPECIFIED  WAVE  DISPERSION  PROPERTIES 

C 

c 

C 

C  THIS  PROGRAM  RUNS  ON  A  PRIME  750  COMPUTER  USING  FORTRAN  77 

C 

C  UNIT  1  IS  THE  KEYBOARD  FOR  AN  INTERACTIVE  TERMINAL 

C 

c 

C 

C  THE  VARIABLES  ENTERED  IN  FREE  FORMAT  ARE: 

C 

C  QO       THE  BASE  FLOW  LEVEL  (CFS) 

C  QP     -  THE  INCREMENTAL  FLOW  ADDED  TO  THE  BASE  FLOW  (CFS) 

C  RQFACT  -  THE  REFERENCE  FLOW  LEVEL,  EXPRESSED  AS  A  DECIMAL 

C  FRACTION  OF  QP 

C  CT       THE  STANDARD  DEVIATION  OF  THE  HYDROGRAPH  GENERATED 

C  (HOURS) 

C  BS     -  IF  EM=5/3,  TOP  WIDTH  OF  THE  CHANNEL  (FEET) 

C  IF  EM=4/3,  SIDE  SLOPE  OF  CHANNEL  (RISE/RUN) 

C  CM     -  MANNINGS  FRICTION  COEFFICIENT 

C  M      -  =0,  FOR  A  RECTANGULAR  CHANNEL 

C  =1,  FOR  A  TRIANGULAR  CHANNEL 

C  SOW    -  CHANNEL  BOTTOM  SLOPE  (FEET/MILE) 

C  DELHR   -  TIME  INTERVAL  BETWEEN  ORDINATES  GENERATED  (HOURS) 

C  LHOUR  -  LAST  HOUR  GENERATED 

C 

c 

C 

C  THE  INPUT  DATA,  INTERMEDIATE  VALUES,  AND  THE  OUTPUT 

C  TIMES  AND  DISCHARGES  ARE  WRITTEN  TO  THE  FILE  DESIGNATED 

C 

c 

C 


CHARACTER*80  FILENAME 
DIMENSION  H(5000) 

PRINT  *, 'ENTER  QO, QP, RQFACT, CT, BS.CM.M, SOM.DELKR, LHOUR' 
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^EAD(1,*)Q0,QP,RQFACT,CT,BS,CM,M,S0M,DELHR,LH0UR 

IF(QO.LT.0.)GO  TO  11 

RSO-SOM/5280. 

IF(M.EQ.l)GO  TO  5 
C 
c 

C 

C  RECTANGULAR  CHANNEL 


EM-5./3. 
RQ-(QO+RQFACT*QP)/BS 

RY-(RQ*CM/(1.486*SQRT(RSO)))**(l./EM) 

RV-RQ/RY 

FN«RV/SQRT(32.2*RY) 

RRF-1.-FN*FN*(1.-EM*(2.-EM)) 

RD«0.5*RQ/RSO*RRF 

GO  TO  8 
C 
c 

C 

C     TRIANGULAR  CHANNEL 

C 

C 

c 


5   EM-4./3. 
Z=BS 
RQ=QO+RQFACT*QP 

RY-((RQ*CW*2.**(2./3.))/(1.486*SQRT(RSO)*Z))**(l./(2.*EM)) 
RV=RQ/(Z*RY*RY) 

FN=RV/SQRT(32.2*RY) 

RRF=1 .-FN*FN*( 1 .-EM*( 2 .-EM) ) 

RD=0 . 5  *RQ / ( RSO* ( 2 . *Z*RY ) ) *RRF 
C 
c 

C 

C     THE  FOLLOWING  CODE  IS  COMMON  TO  BOTH  CHANNEL  SHAPES 

C 

C 

C  ~  """""    ~~ 

8  CTS=CT*360O. 
RRC=EM*RV 

EXO=(((CTS*RRC)**2)*RRC)/(2.*RD) 

CC2=4.*RD 

DXO=EXO*RSO/RY 

RTP-((3.*RD)/RRC**2)*(SQRT(l.+(RRC*EXO/(3.*RD))**2)-l.) 

RHP=EXO/SQRT(4.*3.1416*RD*RTP**3)*EXP(-(EXO-RRC*RTP)**2/(CC2*RTP)) 
RQP=QP/RHP 

CCl=RQP*EXO/SQRT(4.*3.1416*RD) 
NSTEPS=LHOUR/DELHR 

IF(NSTEPS.GT.50O0)NSTEPS»5OO0 
C 
c 
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C 

C     GENERATE  INVERSE  GAUSSIAN  SHAPE  HYDROGRAPH,  H 

C 

C 

c 

TIMINODELHR*3600. 
DO   10   I-l.NSTEPS 
TIME-I*TIMINC 

H(I)-QO+(CCl/(TIME**1.5))*EXP(-(EXO-RRC*TIME)**2/(CC2*TIME)) 
10  CONTINUE 


C 

c 

C     STORE  RESULTS  IN  AN  OUTPUT  FILE 

C 

C 


PRINT  *, 'ENTER  FILE  NAME  FOR  DISCHARGES  AND  TIMES' 
C 

READ (1,1) FILENAME 
1  FORMAT (A80) 

OPEN (FILE-FILENAME , UNIT-1 05 ) 
WRITE(105,500)QO,QP,RQFACT,CT,BS,CM,EM,SOM,DELHR,LHOUR, 

1  RQ,SOM,RSO,RY,RV,FN,RRF,RD,CTS,RRC, 

2  EXO,CC2,DXO,RTP,RHP,RQP,CCl 
C 

500  FORMAT ('QO-', Gil. 4,',  QP-'G11.4,',  RQFACT-' ,G11 .4/ 
1  'CT-',G11.4,',  BS-',G11.4/ 

1  'CM=',G11.4,',  EM-', Gil. 4,',  S0M-',G11.4/ 
1  'DELHR-',G11.4,',  LHOUR-'.Ill/ 
1  'RQ-',G11.4,',  S0M-',G11.4/ 
1  'RS0=',G11.4,',  RY-',G11.4/ 
1  'RV=',G11.4,',  FN-', Gil. 4/ 
1  'RRF-',G11.4,',  RD-',G11.4/ 
1  'CTS=',G11.4,',  RRC=',G11.4/ 
1  'EX0=',G11.4,',  CC2-',G11.4/ 
1  'DX0=',G11.4,',  RTP-',G11.4/ 
1  'RHP-', Gil .4,',  RQP=',G11.4/ 
1  'CC1=',G11.4) 
C 

WRITE(105,501)(H(I),I-1,NSTEPS) 

501  FORKAT(8F10.2) 
C 

11  STOP 
END 


APPENDIX  C 
A  FORTRAN  PROGRAM  FOR  THE  NONLINEAR  STATE  SPACE  ROUTING  MODEL 


C 

c 

C  THIS  PROGRAM  SIMULATES  WITH  THE  NONLINEAR  STATE  SPACE  MODEL 

C 

c 

c 

C  THIS  PROGRAM  RUNS  ON  A  PRIME  750  COMPUTER  UNDER  FORTRAN  77 

C 

C  UNIT  1  IS  THE  KEYBOARD  FOR  AN  INTERACTIVE  TERMINAL 

C 

c 

c 

C  THE  INPUT  READ  FROM  A  FILE  IN  FREE  FORMAT  IS: 

C 

C  AO       THE  PURE  DELAY  AT  THE  REFERENCE  FLOW  LEVEL 

C  NO       THE  NUMBER  OF  RESERVOIRS  AT  THE  REFERENCE  FLOW  LEVEL 

C  (THIS  MAY  BE  A  NONINTEGER  VALUE) 

C  KO     -  THE  PARAMETER  K  AT  THE  REFERENCE  FLOW  LEVEL 

C  QO       THE  REFERENCE  FLOW  LEVEL 

C  M      -  THE  KINEMATIC  WAVE  PARAMETER 

C  DT       THE  TIME  STEP 

C  NSTI   -  THE  NUMBER  OF  VALUES  OF  INFLOW  ENTERED 

C  (DT*NSTI  IS  THE  PERIOD  SIMULATED) 

C  RHO    -  THE  TIME  WEIGHTING  FACTOR  (RANGE  FROM  0  TO  1) 

C  DEBUG  -  =0,  NO  DEBUG  PRINTOUT, 

C  =1,  PRINTOUT  FOR  EACH  RESERVOIR  FOR  EACH  TIME  STEP 

C  INFLOW  -  THE  INFLOW  TIME  SERIES  (NSTI  VALUES) 

C 

c 

C 

C  THE  OUTFLOW  FROM  THE  DOWNSTREAM  RESERVOIR  AND  THE  CORRESPONDING 

C  TIMES  ARE  STORED  IN  ARRAYS  RESOUT  AND  TIMOUT,  RESPECTIVELY. 


CHARACTER*80  FILENAME 
LOGICAL  LINEAR 
INTEGER*4  TOTRES, DEBUG 

REAL*4  K , M , KO , NO .KOLOG , LGOLOG , INFLOW , INFBAR , K I , INFLOG , NO LOG 
C0MM0H/X1X/  F(100,100),G(100),X(100),XP(100),GU(100),K(100), 
1   TLAG(100),AX(100),QLIN(100) 
COMMON/X2X/INFLOW(5001 ) ,RESOUT(5001 ) ,TIMOUT(5001 ) 


NRD=100 
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HWR-101 
C 

PRINT*, 'ENTER  INPUT  FILE  NAME' 
C 

READ (1,1) FILENAME 
1  FORMAT (A80) 

OPEN( ONIT-NRD , FILE-FILENAME ) 
C 

PRINT*, 'ENTER  OUTPUT  FILE  NAME' 
C 

READ(1,1)FILENAME 

OPEN(UNIT-NWR, FILE-FILENAME) 
C 

READ(NRD,*)AO, NO, K0,Q0,M,DT,NSTI,RHO, DEBUG 

IF(N0.EQ.0.)GO  TO  999 
C 

READ(NRD,*)(INFL0W(I),I-1,NSTI) 
C 

NMAX-100 

LINEAR-. FALSE. 

IF(M.EQ.1.)LINEAR-.TRUE. 

IF(LINEAR)GO  TO  400 

TOTLGO-AO+NO/KO 

BLAG-(1.-M)/M 
C 

C     P  =  2/3  FOR  MANNINGS  FRICTION  FORMULA 
C     R  IS  A  FUNCTION  OF  M  FROM  THE  PRISMATIC  SHAPE  FORMULA 
C 

P-2./3. 

R-(l.+P-M)/P 

BN=(R-1.)/M 

BK-(M-2.+R)/M 

QOLOG-LOGIO(QO) 

KOLOG=LOG10(KO) 

AOLOG=LOG10(AO) 

N0LOG=LOG10(N0) 

LGOLOG=LOG10(TOTLGO) 
C 

INFLOG=LOG10 ( INFLOW( 1 ) ) 

TOTLAG=10.**(BLAG*(INFLOG-QOLOG)+LGOLOG) 

N=10.**(BN*(INFLOG-Q0LOG)+N0LOG)  +  0.5 

IF ( N . GT . NMAX ) N-NMAX 

IF(N.LT.1)N=1 

KI=10.**(BK*(INFLOG-Q0LOG)+K0LOG) 

AI=TOTLAG-NI/KI 

GO  TO  401 
C 

400  KI=K0 
N=N0+0.5 
A=A0 

C 

401  S0=INFLOW(l) 
INFBAR=INFL0W(1) 
V0LIN=0. 
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VOLOUTLR-0. 
TOTRES-0 
C 

DO  2  1-1,11 
K(I)-KI 
2  X(I)-S0 

C 

DO  10  J-l.NSTI 

C 

IF(J.GT.l)INFBAR«(INFLOW(J)+INFLOW(J-l))/2. 

C 

IF (LINEAR) GO  TO  405 

C 

C     COMPUTE  AVERAGE  OUTFLOW  FROM  ALL  RESERVOIRS 

C     AS  FLOW  TO  DETERMINE  NEW  NUMBER  OF  RESERVOIRS 

C 

QOBAR-INFBAR 

DO  705  1-1,1 
705  QOBAR"QOBAR+X(I) 

QOBAR-QOBAR/(N+l) 
C 

NOLD-N 
C 

C     COMPUTE  NEW  NUMBER  OF  RESERVOIRS  BASED  ON  QOBAR 
C 

N-10.**(BN*(LOG10(QOBAR)-QOLOG)+NOLOG)  +  0.5 

IF(N.GT.NMAX)N«NMAX 

IF(N.LT.1)N-1 
C 

TOTRES=TOTRES+N 

IF(N.EQ.NOLD)GO  TO  730 
C 

C      IF  NUMBER  OF  RESERVOIRS  HAS  CHANGED  - 
C      INTERPOLATE  OLD  STATES  TO  GET  DISCHARGES  FOR  NEW 
C     NUMBER  OF  RESERVOIRS 
C 

JJ=0 

LCM=N*NOLD 
C 

DO  710  I«1,LCM 

IF((l/NOLD)*NOLD.LT.I)GO  TO  710 

I1=(I-1)/N  +  1 

I2=(I-1)/N0LD  +  1 

QLEFT=INFBAR 

IF(I1 .GT. 1 )QLEFT=X(I1-1 ) 

QRIGHT=X(I1) 

QDIFF=QRIGHT-QLEFT 

XLEFT=N*(I1-1) 

XRIGHT=N*I1 

XWANT=NOLD*I2 

JJ=JJ+1 

QLIN(JJ)=QLEFT+QDIFF*(XWANT-XLEFT)/(XRIGHT-XLEFT) 
710  CONTINUE 
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C     FILL  X  ARRAY  WITH  OUTFLOWS  FROM  EACH  NEW  RESERVOIR 
C 

DO  720  I-l.N 
720  X(I)-QLINU) 
C 

730  ATOT-0. 

DO  3  1-1,1 
C 

IF(I.EQ.l)QBAR-(lNFBAR+X(l))/2. 

IF(I.GT.l)QBAR«(X(I-l)+X(I))/2. 
XXX-LOG10 (QBAR)-QOLOG 
AX ( I ) -1 0 . ** ( BLAG*XXX+A0LOG ) 
K ( I ) -1 0 .** ( BK*XXX+K0L0G ) 
TLAG ( I ) -1 0 . ** ( BLAG*XXX+LG0LOG ) 
ATOT-ATOT+AX(I) 
3  CONTINUE 
A-ATOT/N 
C 

405  IF(.NOT.LINEAR.OR.J.EQ.l)CALL  FILLFG(F,G,NMAX,N,K,DT,RHO,NWR) 

C 

CALL  LDMULT(F,X,XP,N,1,NMAX,1) 


DO  5  I-l.N 
5  GU(I)-G(I)*INFBAR 

CALL  ADD(XP,GU,X,N,1,NMAX,1) 

RESOUT(J)»X(N) 

NDELAY=A/DT+0.5 

ITMOUT=J+NDELAY 
TIMOUT(J)=ITMOUT*DT 

IF(DEBUG.EQ.0)GO  TO  410 


WRITE(NWR,653)J,A,Q0BAR,N,INFBAR,RES0UT(J),TIM0UT(J) 

653  FORMATC  ><><><><><><  TIME  STEP-', 15,',  A«',F10.3,',  QOBAR-', 

1  F10.4,',  N=',I4,',  INFBAR«',F10.3,',  RESOUT(J)-', 

2  F10.4,',  TIMOUT(J)-',F8.2) 
WRITE(NWR,654)(X(I) ,1-1 ,N) 

654  FORMAT ('  DISCHARG-' ,10(1X,F10 .4)/( 10X,10(1X,F10.4) ) ) 
IF(LINEAR)GO  TO  410 
WRITE(KWR,655)(K(I),I-1,N) 

655  FORMATC        K-',10(1X,F10.4)/(10X,10(1X,F10.4))) 
WRITE(NVR,648)(TLAG(I),I=1,N) 

648  FORMATC     TLAG=' ,10( 1X,F10.4)/(10X, 10( 1X.F10.4) )) 
WRITE(KWR,673)(AX(I),I=1 ,N) 

673  FORMATC       AX='  ,10(  1X.F10  .4)/(  10X,  10(  1X.F10.4) ) ) 
IF(N.NE.N0LD)WR1TE(NWR,674)(QLIN(I),I-1,N) 

674  FORMATC     QLIN-' ,10( lX,F10.4)/( 10X, 10( 1X.F10.4))) 
C 

410  VOLIN=VOLIN+INFLOW(J)*DT 

VOLOUTLR=VOL0UTLR+RESOUT(J)*DT 
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C 

10  CONTINUE 
C 

AVGRES-FLOAT(TOTRES)/NSTI 
WRITE(NWR,657)VOLIN,VOLOUTLR,AVGRES 
657  FORMATC  THE  TOTAL  INFLOW  -'.G17.6/ 

2  '  THE  TOTAL  OUTFLOW  FROM  RESERVOIRS  -'.G17.6/ 

3  '  THE  AVERAGE  NUMBER  OF  RESERVOIRS  »  '.F8.3) 
C 

999  STOP 
END 
c 

c 

SUBROUTINE  FILLFG(F,G,MAX,N,K,DT,RHO,NWR) 
C 
C 

c 

C     FILL  F  AND  G  MATRICES  FOR  X(T+1)  -  F*X(T)  +  G*U  EQUATION 

C     FOR  STATES  OF  DISCHARGE 

C 


REAL*4  K 

DIMENSION  F(MAX,MAX),G(MAX),K(MAX) 
C 

C     FILL  MATRIX  F(l,J)  -  I-ROWS, J-COLS 
C 

GAMMA-1 .-RHO 
C 

C     FILL  FROM  1ST  TO  NTH  ROWS 
C     SINCE  K  IS  CONSTANT  FOR  A  ROW 
C 

DO   10    I-l.N 

SQUIGL=1 ./K(I)+DT*RHO 

IF(I.EQ.N)GO  TO  6 

L=I+1 
C 

C  ZERO   TO   THE  RIGHT   OF  MAIN   DIAGONAL 

C 

DO  5   J=L,N 

5  F(I,J)=0. 
C 

C     FILL  IN  VALUE  ON  MAIN  DIAGONAL 
C 

6  F(I,I)=(1./K(I)-DT*GAMMA)/SQUIGL 
IF(I.EQ.l)GO  TO  9 

C 

C     FILL  IN  LEFT  OFF  DIAGONAL 

C 

L=I-1 

F(I,L)=(F(L,L)*RHO+GAMMA)*DT/SQUIGL 

IF(L.EQ.l)GO  TO  9 
C 
C     FILL  TO  THE  LEFT  OF  LEFT  OFF  DIAGONAL 
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C     FILL  FROM  RIGHT  TO  LEFT 
C 

DO  8  JJ-2.L 

J-I-JJ 

F(l,J)-F(L,J)*DT*RHO/SQUIGL 

8  CONTINUE 
C 

C     FILL  IN  G  VECTOR  VALUE 
C 

9  IF(I.EQ.1)G(I)-DT/SQUIGL 
IF(I.GT.1)G(I)-G(L)*DT*RH0/SQUIGL 

10  CONTINUE 
C 

RETURN 
END 


SUBROUTINE  LDMULT(A,B,C,L,M,LDIM,MDIM) 
C 

C     MULTIPLY  TWO  LOWER  DIAGIONAL  MATRICES 
C 


DIMENSION  A(LDIM,LDIM),B(LDIM,MDIM),C(LDIM,MDIM) 

DO  20  I-l.L 

DO  20  J-l.M 

TEMP-0 . 

DO  10  K-1,1 
10  TEMP=TEMP+A(I,K)*B(K,J) 
20  C(I,J)=TEMP 


RETURN 
END 


c 

C 

SUBROUTINE  ADD(A, B,C,N,M,NDIM,MDIM) 
C 

C     ADD  TWO  MATRICES 
C 
c 

C 


DIMENSION  A(NDIM.MDIM) ,B(NDIM,MDIM) ,C(NDIM,MDIM) 
DO  10  1=1, N 
DO  10  J=1,M 
10  C(I,J)=A(I,J)+B(I,J) 
C 

RETURN 
END 
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